#############################################################################################
########### calculate exposure inex for each species ########################################
#############################################################################################

#load packages:
library(R.utils)
library(raster)
library(sp)
library(sf)
library(stringr)
library(circular)

#Main working dir:
MyDir<-paste("/home/jc217070/Maxent",sep="")  #define working directory on HPC
OccDir<-file.path(MyDir,"Maxent_SWDs/spp_swds")
MUOutDir<-paste(MyDir, "/Outputs2_Crossval/permutation.importance", sep="")
IndSppOutDir<-paste(MyDir, "/Outputs3_Final/Final_Individual_Spp", sep="")
PlotDir<-paste(MyDir,"/Maxent_Data/Plots/Plots_WHOSnakes_2020_09_01",sep="")

#set project name:
projectName<-"WHOSnakes"

#define 'for' variables:
crit<-c("permutation.importance")   #c("permutation.importance","contribution","Test.gain.without","Test.gain.with.only","AUC.without","AUC.with.only")
timecats<-c("Current","2050_Median","2050_Q90","2050_Q10","2090_Median","2090_Q90","2090_Q10")
#timecats<-c("Current","2050_Median","2090_Median")
#timecats<-c("Current")

#read species reference list:
DatSum<-read.table(file=file.path(paste(MyDir, "/Maxent_LUTables/",projectName,".csv",sep="")), header=TRUE, sep=",") #load data summary for all species
DatSum<-DatSum[DatSum$Listing_comment!="not listed",]
spp<-DatSum$gen_sp_subtax[DatSum$Listing_comment!="not listed"]
sppA<-DatSum$Accepted_Species[DatSum$Listing_comment!="not listed"]
#spSTART<-which (spp=="Micrurus_mertensi")







#DatSumMod<-read.table(file=file.path(paste(MyDir, "/Maxent_LUTables/",projectName,"_SupMatFin.csv",sep="")), header=TRUE, sep=",") #load data summary for all species
DatSumMod<-DatSum[,c("SSN","Family","Accepted_Species","Listing_Name","ModUnit","MergeSpecies","MergeShort","country")]
DatSumMod<-DatSumMod[DatSumMod$Accepted_Species %in% sppA,]


# make extra columns to fill
DatSumMod$SPrecs <- NA
DatSumMod$MUrecs <- NA
DatSumMod$TrainingAUC <- NA
DatSumMod$TestAUC <- NA

DatSumMod$RangesizeSQKM_Current <- NA
DatSumMod$SuitabilitySum_Current <- NA

DatSumMod$RangesizeSQKM_2050_Median <- NA
DatSumMod$SuitabilitySum_2050_Median <- NA
DatSumMod$RangesizeSQKM_2090_Median <- NA
DatSumMod$SuitabilitySum_2090_Median <- NA

DatSumMod$RangesizeSQKM_2050_Q10 <- NA
DatSumMod$SuitabilitySum_2050_Q10 <- NA
DatSumMod$RangesizeSQKM_2090_Q10 <- NA
DatSumMod$SuitabilitySum_2090_Q10 <- NA

DatSumMod$RangesizeSQKM_2050_Q90 <- NA
DatSumMod$SuitabilitySum_2050_Q90 <- NA
DatSumMod$RangesizeSQKM_2090_Q90 <- NA
DatSumMod$SuitabilitySum_2090_Q90 <- NA

DatSumMod$lat_Current <- NA
DatSumMod$lat_2050_Median <- NA
DatSumMod$lat_2090_Median <- NA
DatSumMod$lat_2050_Q10 <- NA
DatSumMod$lat_2050_Q90 <- NA
DatSumMod$lat_2090_Q10 <- NA
DatSumMod$lat_2090_Q90 <- NA

DatSumMod$long_Current <- NA
DatSumMod$long_2050_Median <- NA
DatSumMod$long_2090_Median <- NA
DatSumMod$long_2050_Q10 <- NA
DatSumMod$long_2050_Q90 <- NA
DatSumMod$long_2090_Q10 <- NA
DatSumMod$long_2090_Q90 <- NA




DatSumMod$Risk_Current <- NA
DatSumMod$Exposure_Current <- NA
DatSumMod$Impact_Current <- NA

DatSumMod$Risk_2050_Median <- NA
DatSumMod$Exposure_2050_Median <- NA
DatSumMod$Impact_2050_Median <- NA
DatSumMod$Risk_2090_Median <- NA
DatSumMod$Exposure_2090_Median <- NA
DatSumMod$Impact_2090_Median <- NA

DatSumMod$Risk_2050_Q10 <- NA
DatSumMod$Exposure_2050_Q10 <- NA
DatSumMod$Impact_2050_Q10 <- NA
DatSumMod$Risk_2090_Q10 <- NA
DatSumMod$Exposure_2090_Q10 <- NA
DatSumMod$Impact_2090_Q10 <- NA

DatSumMod$Risk_2050_Q90 <- NA
DatSumMod$Exposure_2050_Q90 <- NA
DatSumMod$Impact_2050_Q90 <- NA
DatSumMod$Risk_2090_Q90 <- NA
DatSumMod$Exposure_2090_Q90 <- NA
DatSumMod$Impact_2090_Q90 <- NA


# read LN people per cell & % agriculture
if(file.exists(paste("/home/jc217070/GIS_and_raw_Layers/final_layers/wpppcln.asc.gz")) &
   !file.exists(paste("/home/jc217070/GIS_and_raw_Layers/final_layers/wpppcln.asc"))){
  gunzip(filename=paste("/home/jc217070/GIS_and_raw_Layers/final_layers/wpppcln.asc.gz"),
         destname=paste("/home/jc217070/GIS_and_raw_Layers/final_layers/wpppcln.asc"),
         remove=FALSE, overwrite=TRUE)
}
if(file.exists(paste("/home/jc217070/GIS_and_raw_Layers/final_layers/percag.asc.gz")) &
   !file.exists(paste("/home/jc217070/GIS_and_raw_Layers/final_layers/percag.asc"))){
  gunzip(filename=paste("/home/jc217070/GIS_and_raw_Layers/final_layers/percag.asc.gz"),
         destname=paste("/home/jc217070/GIS_and_raw_Layers/final_layers/percag.asc"),
         remove=FALSE, overwrite=TRUE)
}

ppc<-raster(paste("/home/jc217070/GIS_and_raw_Layers/final_layers/wpppcln.asc"), crs="+init=epsg:4326")
pag<-raster(paste("/home/jc217070/GIS_and_raw_Layers/final_layers/percag.asc"), crs="+init=epsg:4326")

worldmap <- st_read(paste(PlotDir,"/worldmap/WHO_countries.shp",sep=""))
WHOpolies <- st_read(paste(MyDir,"/WHOSnakesShapes/SSN_CTY_CAT_BRG_FEA_IN_Mar2023.shp",sep=""))


for (sp in spp){
#for (sp in spp[spSTART:length(spp)]){

  print(paste(Sys.time(), "appending table for", sp, sep=" "))

  print(paste(Sys.time(), "model stats", sep=" "))
  
  n<-which(spp==sp) #gets row number for current species
  mapsp<-DatSum$ModUnit[n]
  
  if(file.exists(paste(OccDir, "/", sp,".csv",sep=""))){
  sprecs<-read.table(file=file.path(paste(OccDir, "/", sp,".csv",sep="")), header=TRUE, sep=",")
  NSPrecs<-length(sprecs[,1])
  DatSumMod$SPrecs [n] <- NSPrecs
  }
  
  maprecs<-read.table(file=file.path(paste(OccDir, "/ModUnit_", gsub(" ","_",mapsp),".csv",sep="")), header=TRUE, sep=",")
  NMUrecs<-length(maprecs[,1])
  DatSumMod$MUrecs [n] <- NMUrecs
  
  results <-read.csv(file=paste(MUOutDir, "/ModUnit_", gsub(" ","_",mapsp),"/maxentResults.csv", sep=""), header=TRUE)
  DatSumMod$TrainingAUC [n] <- results$Training.AUC[length(results$Training.AUC)]
  DatSumMod$TestAUC [n] <- results$Test.AUC[length(results$Test.AUC)]

# calculate other ?!!!stats!!!?, 

for (cat in timecats){
  
  print(paste(Sys.time(), "SDM stats", sep=" "))
  
  outfile<-paste(IndSppOutDir,"/",sp,"_",cat, "_CropThreshCDCut.asc",sep="")
  
  # unzip raster
  if(file.exists(paste(outfile,".gz", sep="")) & !file.exists(outfile)){
  gunzip(filename=paste(outfile,".gz", sep=""),
         destname=paste(outfile),
         remove=FALSE, overwrite=TRUE)
  }
  
  # read raster
  SDM<-raster(paste(outfile), crs="+init=epsg:4326")
  
  
  # get all values
  myvals <- values(SDM)
  myvals<- myvals[myvals>0]
  myvals<- myvals[!is.na(myvals)]
  mysum<-sum(myvals)
  myrange <- length(myvals)
  
  
  mycoords<-as.data.frame(xyFromCell(SDM, which(SDM[]>0)))
  myX<-mean(mycoords$x)
  myY<-mean(mycoords$y)
  
  # append table with sum (total suitability) & range size
  mycolname<-paste("RangesizeSQKM_",cat,sep="")
  myc<-which(names(DatSumMod)==mycolname)
  DatSumMod [n,myc] <- myrange
  
  mycolname<-paste("SuitabilitySum_",cat,sep="")
  myc<-which(names(DatSumMod)==mycolname)
  DatSumMod [n,myc] <- mysum
  
  mycolname<-paste("lat_",cat,sep="")
  myc<-which(names(DatSumMod)==mycolname)
  DatSumMod [n,myc] <- myY
  
  mycolname<-paste("long_",cat,sep="")
  myc<-which(names(DatSumMod)==mycolname)
  DatSumMod [n,myc] <- myX
  
  ##############################################
  
  print(paste(Sys.time(), "Risk Index", sep=" "))
  
#Risk Index
Risk<-pag*SDM

#SAVE!!!!#
writeRaster(Risk, filename=paste(IndSppOutDir,"/",sp,"_",cat,"_Risk.asc",sep=""), format="ascii", overwrite=TRUE)
gzip(filename=paste(IndSppOutDir,"/",sp,"_",cat,"_Risk.asc",sep=""),overwrite=TRUE, remove=TRUE)

# get all values & sum
myvals <- values(Risk)
myvals<- myvals[myvals>0]
myvals<- myvals[!is.na(myvals)]
mysum<-sum(myvals)
# append table row with sum (risk index)
mycolname<-paste("Risk_",cat,sep="")
myc<-which(names(DatSumMod)==mycolname)
DatSumMod [n,myc] <- mysum

print(paste(Sys.time(), "Exposure Index", sep=" "))

#Exposure Index
Exp<-ppc*SDM


#SAVE!!!!#
writeRaster(Exp, filename=paste(IndSppOutDir,"/",sp,"_",cat,"_Exp.asc",sep=""), format="ascii", overwrite=TRUE)
gzip(filename=paste(IndSppOutDir,"/",sp,"_",cat,"_Exp.asc",sep=""),overwrite=TRUE, remove=TRUE)


myvals <- values(Exp)
myvals<- myvals[myvals>0]
myvals<- myvals[!is.na(myvals)]
mysum<-sum(myvals)
DatSumMod [n,myc+1] <- mysum

print(paste(Sys.time(), "Impact Index", sep=" "))

#Impact Index
Imp<-Exp*pag


#SAVE!!!!#
writeRaster(Imp, filename=paste(IndSppOutDir,"/",sp,"_",cat,"_Imp.asc",sep=""), format="ascii", overwrite=TRUE)
gzip(filename=paste(IndSppOutDir,"/",sp,"_",cat,"_Imp.asc",sep=""),overwrite=TRUE, remove=TRUE)


myvals <- values(Imp)
myvals<- myvals[myvals>0]
myvals<- myvals[!is.na(myvals)]
mysum<-sum(myvals)
DatSumMod [n,myc+2] <- mysum



#remove unzipped raster or zip
if(file.exists(paste(outfile,".gz", sep="")) & file.exists(outfile)){
  file.remove(outfile)
}
if(!file.exists(paste(outfile,".gz", sep="")) & file.exists(outfile)){
  gzip(filename=paste(outfile))
}

#################################
#################################
#read in range loss/gain raster:

outfile<-paste(IndSppOutDir,"/",sp,"_",cat, "_LossVsGain.asc",sep="")

# unzip raster
if(file.exists(paste(outfile,".gz", sep="")) & !file.exists(outfile)){
  gunzip(filename=paste(outfile,".gz", sep=""),
         destname=paste(outfile),
         remove=FALSE, overwrite=TRUE)
}

# read raster
SDM<-raster(paste(outfile), crs="+init=epsg:4326")

# get all values
myvals <- values(SDM)

#calculate gained range 
myvals<- myvals[myvals>0]
myvals<- myvals[!is.na(myvals)]
mygain <- length(myvals)

# get all values
myvals <- values(SDM)

#calculate lost range
myvals<- myvals[myvals<0]
myvals<- myvals[!is.na(myvals)]
myloss <- length(myvals)


#add to summary file
mycolname<-paste("Gain_",cat,sep="")
myc<-which(names(DatSumMod)==mycolname)
DatSumMod [n,myc] <- mygain

mycolname<-paste("Loss_",cat,sep="")
myc<-which(names(DatSumMod)==mycolname)
DatSumMod [n,myc] <- myloss

#remove unzipped raster or zip
if(file.exists(paste(outfile,".gz", sep="")) & file.exists(outfile)){
  file.remove(outfile)
}
if(!file.exists(paste(outfile,".gz", sep="")) & file.exists(outfile)){
  gzip(filename=paste(outfile))
}

#################################
#################################

#write copy of updated table
write.table(x=DatSumMod, file=paste(MyDir, "/Maxent_LUTables/",projectName,"_SupMat.csv",sep=""),
            na = "NA", row.names = FALSE, col.names = TRUE, sep=",")
}}



DatSumMod$RangesizeSQKM_CHANGE_2050_Median  <-  DatSumMod$RangesizeSQKM_2050_Median - DatSumMod$RangesizeSQKM_Current
DatSumMod$RangesizeSQKM_CHANGE_2050_Q10     <-  DatSumMod$RangesizeSQKM_2050_Q10 - DatSumMod$RangesizeSQKM_Current
DatSumMod$RangesizeSQKM_CHANGE_2050_Q90     <-  DatSumMod$RangesizeSQKM_2050_Q90 - DatSumMod$RangesizeSQKM_Current
DatSumMod$RangesizeSQKM_CHANGE_2090_Median  <-  DatSumMod$RangesizeSQKM_2090_Median - DatSumMod$RangesizeSQKM_Current
DatSumMod$RangesizeSQKM_CHANGE_2090_Q10     <-  DatSumMod$RangesizeSQKM_2090_Q10 - DatSumMod$RangesizeSQKM_Current
DatSumMod$RangesizeSQKM_CHANGE_2090_Q90     <- DatSumMod$RangesizeSQKM_2090_Q90 - DatSumMod$RangesizeSQKM_Current

DatSumMod$SuitabilitySum_CHANGE_2050_Median <- DatSumMod$SuitabilitySum_2050_Median - DatSumMod$SuitabilitySum_Current
DatSumMod$SuitabilitySum_CHANGE_2050_Q10    <- DatSumMod$SuitabilitySum_2050_Q10 - DatSumMod$SuitabilitySum_Current
DatSumMod$SuitabilitySum_CHANGE_2050_Q90    <- DatSumMod$SuitabilitySum_2050_Q90 - DatSumMod$SuitabilitySum_Current
DatSumMod$SuitabilitySum_CHANGE_2090_Median <- DatSumMod$SuitabilitySum_2090_Median - DatSumMod$RangesizeSQKM_Current
DatSumMod$SuitabilitySum_CHANGE_2090_Q10    <- DatSumMod$SuitabilitySum_2090_Q10 - DatSumMod$RangesizeSQKM_Current
DatSumMod$SuitabilitySum_CHANGE_2090_Q90    <- DatSumMod$SuitabilitySum_2090_Q90 - DatSumMod$RangesizeSQKM_Current

DatSumMod$Risk_CHANGE_2050_Median           <-  DatSumMod$Risk_2050_Median -  DatSumMod$Risk_Current
DatSumMod$Risk_CHANGE_2050_Q10              <-  DatSumMod$Risk_2050_Q10 -  DatSumMod$Risk_Current
DatSumMod$Risk_CHANGE_2050_Q90              <-  DatSumMod$Risk_2050_Q90 -  DatSumMod$Risk_Current
DatSumMod$Risk_CHANGE_2090_Median           <- DatSumMod$Risk_2090_Median -  DatSumMod$Risk_Current
DatSumMod$Risk_CHANGE_2090_Q10              <- DatSumMod$Risk_2090_Q10 -  DatSumMod$Risk_Current
DatSumMod$Risk_CHANGE_2090_Q90              <- DatSumMod$Risk_2090_Q90 -  DatSumMod$Risk_Current

DatSumMod$Exposure_CHANGE_2050_Median       <- DatSumMod$Exposure_2050_Median -  DatSumMod$Exposure_Current
DatSumMod$Exposure_CHANGE_2050_Q10          <- DatSumMod$Exposure_2050_Q10 -  DatSumMod$Exposure_Current
DatSumMod$Exposure_CHANGE_2050_Q90          <- DatSumMod$Exposure_2050_Q90 -  DatSumMod$Exposure_Current
DatSumMod$Exposure_CHANGE_2090_Median       <- DatSumMod$Exposure_2090_Median -  DatSumMod$Exposure_Current
DatSumMod$Exposure_CHANGE_2090_Q10          <- DatSumMod$Exposure_2090_Q10 -  DatSumMod$Exposure_Current
DatSumMod$Exposure_CHANGE_2090_Q90          <- DatSumMod$Exposure_2090_Q90 -  DatSumMod$Exposure_Current

DatSumMod$Impact_CHANGE_2050_Median         <- DatSumMod$Impact_2050_Median -  DatSumMod$Impact_Current
DatSumMod$Impact_CHANGE_2050_Q10            <- DatSumMod$Impact_2050_Q10 -  DatSumMod$Impact_Current
DatSumMod$Impact_CHANGE_2050_Q90            <- DatSumMod$Impact_2050_Q90 -  DatSumMod$Impact_Current
DatSumMod$Impact_CHANGE_2090_Median         <- DatSumMod$Impact_2090_Median -  DatSumMod$Impact_Current
DatSumMod$Impact_CHANGE_2090_Q10            <- DatSumMod$Impact_2090_Q10 -  DatSumMod$Impact_Current
DatSumMod$Impact_CHANGE_2090_Q90            <- DatSumMod$Impact_2090_Q90 -  DatSumMod$Impact_Current

#calculate change vectors (distance and direction)
#set up variables:
R<- 6371e3
latpiCurrent = DatSumMod$lat_Current * pi/180
latpi2050Med = DatSumMod$lat_2050_Median * pi/180
latpi2090Med = DatSumMod$lat_2090_Median * pi/180
latpi2050Q10 = DatSumMod$lat_2050_Q10 * pi/180
latpi2050Q90 = DatSumMod$lat_2050_Q90 * pi/180
latpi2090Q10 = DatSumMod$lat_2090_Q10 * pi/180
latpi2090Q90 = DatSumMod$lat_2090_Q90 * pi/180
longpiCurrent = DatSumMod$long_Current * pi/180
longpi2050Med = DatSumMod$long_2050_Median * pi/180
longpi2090Med = DatSumMod$long_2090_Median * pi/180
longpi2050Q10 = DatSumMod$long_2050_Q10 * pi/180
longpi2050Q90 = DatSumMod$long_2050_Q90 * pi/180
longpi2090Q10 = DatSumMod$long_2090_Q10 * pi/180
longpi2090Q90 = DatSumMod$lat_2090_Q90 * pi/180

#calculate distances:
DatSumMod$DistM_2050_Median <-  R * (2*atan2(sqrt     
                                             (sin((latpi2050Med-latpiCurrent)/2) * sin((latpi2050Med-latpiCurrent)/2) +  cos(latpiCurrent) * cos(latpi2050Med) * sin((longpi2050Med-longpiCurrent)/2) * sin((longpi2050Med-longpiCurrent)/2))   
                                             , sqrt(1-  
                                                      (sin((latpi2050Med-latpiCurrent)/2) * sin((latpi2050Med-latpiCurrent)/2) +  cos(latpiCurrent) * cos(latpi2050Med) *  sin((longpi2050Med-longpiCurrent)/2) * sin((longpi2050Med-longpiCurrent)/2)) 
                                             )))

DatSumMod$DistM_2050_Q10 <-  R * (2*atan2(sqrt     
                                             (sin((latpi2050Q10-latpiCurrent)/2) * sin((latpi2050Q10-latpiCurrent)/2) +  cos(latpiCurrent) * cos(latpi2050Q10) * sin((longpi2050Q10-longpiCurrent)/2) * sin((longpi2050Q10-longpiCurrent)/2))   
                                             , sqrt(1-  
                                                      (sin((latpi2050Q10-latpiCurrent)/2) * sin((latpi2050Q10-latpiCurrent)/2) +  cos(latpiCurrent) * cos(latpi2050Q10) *  sin((longpi2050Q10-longpiCurrent)/2) * sin((longpi2050Q10-longpiCurrent)/2)) 
                                             )))

DatSumMod$DistM_2050_Q90 <-  R * (2*atan2(sqrt     
                                             (sin((latpi2050Q90-latpiCurrent)/2) * sin((latpi2050Q90-latpiCurrent)/2) +  cos(latpiCurrent) * cos(latpi2050Q90) * sin((longpi2050Q90-longpiCurrent)/2) * sin((longpi2050Q90-longpiCurrent)/2))   
                                             , sqrt(1-  
                                                      (sin((latpi2050Q90-latpiCurrent)/2) * sin((latpi2050Q90-latpiCurrent)/2) +  cos(latpiCurrent) * cos(latpi2050Q90) *  sin((longpi2050Q90-longpiCurrent)/2) * sin((longpi2050Q90-longpiCurrent)/2)) 
                                             )))

DatSumMod$DistM_2090_Median <-  R * (2*atan2(sqrt     
                                             (sin((latpi2090Med-latpiCurrent)/2) * sin((latpi2090Med-latpiCurrent)/2) +  cos(latpiCurrent) * cos(latpi2090Med) * sin((longpi2090Med-longpiCurrent)/2) * sin((longpi2090Med-longpiCurrent)/2))   
                                             , sqrt(1-  
                                                      (sin((latpi2090Med-latpiCurrent)/2) * sin((latpi2090Med-latpiCurrent)/2) +  cos(latpiCurrent) * cos(latpi2090Med) *  sin((longpi2090Med-longpiCurrent)/2) * sin((longpi2090Med-longpiCurrent)/2)) 
                                             )))

DatSumMod$DistM_2090_Q10 <-  R * (2*atan2(sqrt     
                                          (sin((latpi2090Q10-latpiCurrent)/2) * sin((latpi2090Q10-latpiCurrent)/2) +  cos(latpiCurrent) * cos(latpi2090Q10) * sin((longpi2090Q10-longpiCurrent)/2) * sin((longpi2090Q10-longpiCurrent)/2))   
                                          , sqrt(1-  
                                                   (sin((latpi2090Q10-latpiCurrent)/2) * sin((latpi2090Q10-latpiCurrent)/2) +  cos(latpiCurrent) * cos(latpi2090Q10) *  sin((longpi2090Q10-longpiCurrent)/2) * sin((longpi2090Q10-longpiCurrent)/2)) 
                                          )))

DatSumMod$DistM_2090_Q90 <-  R * (2*atan2(sqrt     
                                          (sin((latpi2090Q90-latpiCurrent)/2) * sin((latpi2090Q90-latpiCurrent)/2) +  cos(latpiCurrent) * cos(latpi2090Q90) * sin((longpi2090Q90-longpiCurrent)/2) * sin((longpi2090Q90-longpiCurrent)/2))   
                                          , sqrt(1-  
                                                   (sin((latpi2090Q90-latpiCurrent)/2) * sin((latpi2090Q90-latpiCurrent)/2) +  cos(latpiCurrent) * cos(latpi2090Q90) *  sin((longpi2090Q90-longpiCurrent)/2) * sin((longpi2090Q90-longpiCurrent)/2)) 
                                          )))

#calculate directions:
DatSumMod$DirDeg_2050_Median = (180/pi) * atan2(sin(longpi2050Med-longpiCurrent) * cos(latpi2050Med) , cos(latpiCurrent)*sin(latpi2050Med) - sin(latpiCurrent)*cos(latpi2050Med)*cos(longpi2050Med-longpiCurrent))
DatSumMod$DirDeg_2050_Q10 = (180/pi) * atan2(sin(longpi2050Q10-longpiCurrent) * cos(latpi2050Q10) , cos(latpiCurrent)*sin(latpi2050Q10) - sin(latpiCurrent)*cos(latpi2050Q10)*cos(longpi2050Q10-longpiCurrent))
DatSumMod$DirDeg_2050_Q90 = (180/pi) * atan2(sin(longpi2050Q90-longpiCurrent) * cos(latpi2050Q90) , cos(latpiCurrent)*sin(latpi2050Q90) - sin(latpiCurrent)*cos(latpi2050Q90)*cos(longpi2050Q90-longpiCurrent))
DatSumMod$DirDeg_2090_Median = (180/pi) * atan2(sin(longpi2090Med-longpiCurrent) * cos(latpi2090Med) , cos(latpiCurrent)*sin(latpi2090Med) - sin(latpiCurrent)*cos(latpi2090Med)*cos(longpi2090Med-longpiCurrent))
DatSumMod$DirDeg_2090_Q10 = (180/pi) * atan2(sin(longpi2090Q10-longpiCurrent) * cos(latpi2090Q10) , cos(latpiCurrent)*sin(latpi2090Q10) - sin(latpiCurrent)*cos(latpi2090Q10)*cos(longpi2090Q10-longpiCurrent))
DatSumMod$DirDeg_2090_Q90 = (180/pi) * atan2(sin(longpi2090Q90-longpiCurrent) * cos(latpi2090Q90) , cos(latpiCurrent)*sin(latpi2090Q90) - sin(latpiCurrent)*cos(latpi2090Q90)*cos(longpi2090Q90-longpiCurrent))
# + = east, - = west
# 0 = north, 90 = east, -90 = west, 180/-180 =south



# calculate summary stats for each value across all species for each time category

DatSumMod<-DatSumMod[DatSumMod$Accepted_Species!="Average",]
DatSumMod[length(DatSumMod[,1])+1,"Accepted_Species"]<- "Average"

for (col in   c(9:(length(DatSumMod)-6))  ){
  
  DatSumMod[length(DatSumMod[,1]),col]  <-mean(DatSumMod[1:(length(DatSumMod$Accepted_Species)-1),col],na.rm=TRUE)
  
}

#do the circular mean for directions...
for (col in   c((length(DatSumMod)-5):length(DatSumMod))){
  #DatSumMod[length(DatSumMod[,1]),col]  <-mean(DatSumMod[,col],na.rm=TRUE)
  
  mydata  <-DatSumMod[1:(length(DatSumMod$Accepted_Species)-1),col]
  mydataC<-circular(mydata, type = "directions", 
                    units = "degrees",
                    modulo = "asis", 
                    rotation = "clock")
  
  DatSumMod[length(DatSumMod[,1]),col]<-mean.circular(mydata,na.rm=TRUE)
}

write.table(x=DatSumMod, file=paste(MyDir, "/Maxent_LUTables/",projectName,"_SupMat.csv",sep=""),
            na = "NA", row.names = FALSE, col.names = TRUE, sep=",")

# write.table(x=DatSumMod, file=paste(MyDir, "/Maxent_LUTables/",projectName,"_SupMatFin.csv",sep=""),
#             na = "NA", row.names = FALSE, col.names = TRUE, sep=",")
# 
# 










#################################################################################
#################################################################################
# plot trends (individual species lines, and overall line)

DatSumMod<-read.table(file=paste(MyDir, "/Maxent_LUTables/",projectName,"_SupMatFinFin.csv",sep=""), header=TRUE, sep=",") #load data summary for all species


###################################
# 0. Individual Species Plots #####
###################################


myX<-c(1985,2050,2090)
#####################################################################
# DatSumMod$RangesizeSQKM_Current[DatSumMod$RangesizeSQKM_Current==0]<-0.01
# DatSumMod$SuitabilitySum_Current[DatSumMod$SuitabilitySum_Current==0]<-0.01
# DatSumMod$Exposure_Current[DatSumMod$Exposure_Current==0]<-0.01
# DatSumMod$Impact_Current[DatSumMod$Impact_Current==0]<-0.01

for (col in c(13:26,41:61)){
  DatSumMod[,col]<-  ifelse(DatSumMod[,col]==0,1,DatSumMod[,col])
}
#####################################################################

for (sp in spp){
  
  print(sp)
  
  n<-which(spp==sp)
  
{  #####################################################################
  # #read lost vs gained files for calculations
  # outfile<-paste(IndSppOutDir,"/",sp,"_","2050_Median_LossVsGain.asc",sep="")
  # # unzip raster
  # if(file.exists(paste(outfile,".gz", sep="")) & !file.exists(outfile)){
  #   gunzip(filename=paste(outfile,".gz", sep=""),
  #          destname=paste(outfile),
  #          remove=FALSE, overwrite=TRUE)
  # }
  # LvG2050<-raster(paste(outfile), crs="+init=epsg:4326")
  # 
  # outfile<-paste(IndSppOutDir,"/",sp,"_","2090_Median_LossVsGain.asc",sep="")
  # # unzip raster
  # if(file.exists(paste(outfile,".gz", sep="")) & !file.exists(outfile)){
  #   gunzip(filename=paste(outfile,".gz", sep=""),
  #          destname=paste(outfile),
  #          remove=FALSE, overwrite=TRUE)
  # }
  # LvG2090<-raster(paste(outfile), crs="+init=epsg:4326")
}
  RangeMed<-c(DatSumMod$RangesizeSQKM_Current[n]/DatSumMod$RangesizeSQKM_Current[n],
              DatSumMod$RangesizeSQKM_2050_Median[n]/DatSumMod$RangesizeSQKM_Current[n],
              DatSumMod$RangesizeSQKM_2090_Median[n]/DatSumMod$RangesizeSQKM_Current[n])
  RLossMed<-c(1,
              1-DatSumMod$Loss_2050_Median[n]/DatSumMod$RangesizeSQKM_Current[n],
              1-DatSumMod$Loss_2090_Median[n]/DatSumMod$RangesizeSQKM_Current[n])
  RGainMed<-c(1,
              1+DatSumMod$Gain_2050_Median[n]/DatSumMod$RangesizeSQKM_Current[n],
              1+DatSumMod$Gain_2090_Median[n]/DatSumMod$RangesizeSQKM_Current[n])
  
  RangeQ10<-c(DatSumMod$RangesizeSQKM_Current[n]/DatSumMod$RangesizeSQKM_Current[n],
              DatSumMod$RangesizeSQKM_2050_Q10[n]/DatSumMod$RangesizeSQKM_Current[n],
              DatSumMod$RangesizeSQKM_2090_Q10[n]/DatSumMod$RangesizeSQKM_Current[n])
  RLossQ10<-c(1,
              1-DatSumMod$Loss_2050_Q10[n]/DatSumMod$RangesizeSQKM_Current[n],
              1-DatSumMod$Loss_2090_Q10[n]/DatSumMod$RangesizeSQKM_Current[n])
  RGainQ10<-c(1,
              1+DatSumMod$Gain_2050_Q10[n]/DatSumMod$RangesizeSQKM_Current[n],
              1+DatSumMod$Gain_2090_Q10[n]/DatSumMod$RangesizeSQKM_Current[n])
  
  RangeQ90<-c(DatSumMod$RangesizeSQKM_Current[n]/DatSumMod$RangesizeSQKM_Current[n],
              DatSumMod$RangesizeSQKM_2050_Q90[n]/DatSumMod$RangesizeSQKM_Current[n],
              DatSumMod$RangesizeSQKM_2090_Q90[n]/DatSumMod$RangesizeSQKM_Current[n])
  RLossQ90<-c(1,
              1-DatSumMod$Loss_2050_Q90[n]/DatSumMod$RangesizeSQKM_Current[n],
              1-DatSumMod$Loss_2090_Q90[n]/DatSumMod$RangesizeSQKM_Current[n])
  RGainQ90<-c(1,
              1+DatSumMod$Gain_2050_Q90[n]/DatSumMod$RangesizeSQKM_Current[n],
              1+DatSumMod$Gain_2090_Q90[n]/DatSumMod$RangesizeSQKM_Current[n])
  #####################################################################
  SuitMed<-c(DatSumMod$SuitabilitySum_Current[n]/DatSumMod$SuitabilitySum_Current[n],
             DatSumMod$SuitabilitySum_2050_Median[n]/DatSumMod$SuitabilitySum_Current[n],
             DatSumMod$SuitabilitySum_2090_Median[n]/DatSumMod$SuitabilitySum_Current[n])
  SuitQ10<-c(DatSumMod$SuitabilitySum_Current[n]/DatSumMod$SuitabilitySum_Current[n],
             DatSumMod$SuitabilitySum_2050_Q10[n]/DatSumMod$SuitabilitySum_Current[n],
             DatSumMod$SuitabilitySum_2090_Q10[n]/DatSumMod$SuitabilitySum_Current[n])
  SuitQ90<-c(DatSumMod$SuitabilitySum_Current[n]/DatSumMod$SuitabilitySum_Current[n],
             DatSumMod$SuitabilitySum_2050_Q90[n]/DatSumMod$SuitabilitySum_Current[n],
             DatSumMod$SuitabilitySum_2090_Q90[n]/DatSumMod$SuitabilitySum_Current[n])
  
{
  
  
  
  # # calculate suitability lost & gained
  # outfile<-paste(IndSppOutDir,"/",sp,"_","Current_CropThreshCDCut.asc",sep="")
  # # unzip raster
  # if(file.exists(paste(outfile,".gz", sep="")) & !file.exists(outfile)){
  #   gunzip(filename=paste(outfile,".gz", sep=""),
  #          destname=paste(outfile),
  #          remove=FALSE, overwrite=TRUE)
  # }
  # Suit1985<-raster(paste(outfile), crs="+init=epsg:4326")
  # 
  # outfile<-paste(IndSppOutDir,"/",sp,"_","2050_Median_CropThreshCDCut.asc",sep="")
  # # unzip raster
  # if(file.exists(paste(outfile,".gz", sep="")) & !file.exists(outfile)){
  #   gunzip(filename=paste(outfile,".gz", sep=""),
  #          destname=paste(outfile),
  #          remove=FALSE, overwrite=TRUE)
  # }
  # Suit2050<-raster(paste(outfile), crs="+init=epsg:4326")
  # 
  # outfile<-paste(IndSppOutDir,"/",sp,"_","2090_Median_CropThreshCDCut.asc",sep="")
  # # unzip raster
  # if(file.exists(paste(outfile,".gz", sep="")) & !file.exists(outfile)){
  #   gunzip(filename=paste(outfile,".gz", sep=""),
  #          destname=paste(outfile),
  #          remove=FALSE, overwrite=TRUE)
  # }
  # Suit2090<-raster(paste(outfile), crs="+init=epsg:4326")
  # 
  # myras<-LvG2050*Suit2050 #(anything>0)
  # myvals <- values(myras)
  # myvals<- myvals[myvals>0]
  # myvals<- myvals[!is.na(myvals)]
  # suitgain2050<-sum(myvals)
  # 
  # myras<-LvG2090*Suit2090 #(anything>0)
  # myvals <- values(myras)
  # myvals<- myvals[myvals>0]
  # myvals<- myvals[!is.na(myvals)]
  # suitgain2090<-sum(myvals)
  # 
  # myras<-LvG2050*Suit1985 #(anything<0)
  # myvals <- values(myras)
  # myvals<- myvals[myvals>0]
  # myvals<- myvals[!is.na(myvals)]
  # suitloss2050<-sum(myvals)
  # 
  # myras<-LvG2090*Suit1985 #(anything<0)
  # myvals <- values(myras)
  # myvals<- myvals[myvals>0]
  # myvals<- myvals[!is.na(myvals)]
  # suitloss2090<-sum(myvals)
  # 
  # SLossMed<-c(0,
  #             suitloss2050/DatSumMod$SuitabilitySum_Current[n],
  #             suitloss2090/DatSumMod$SuitabilitySum_Current[n])
  # SLossMed[SLossMed==0]<-0.01 #avoid dividing by 0 for % change graphs
  # 
  # if(file.exists(paste(IndSppOutDir,"/",sp,"_","2050_Median_CropThreshCDCut.asc.gz",sep="")) & 
  #    !file.exists(paste(IndSppOutDir,"/",sp,"_","2050_Median_CropThreshCDCut.asc",sep=""))){
  #   file.remove(paste(IndSppOutDir,"/",sp,"_","2050_Median_CropThreshCDCut.asc",sep=""))
  # }
  # if(file.exists(paste(IndSppOutDir,"/",sp,"_","2090_Median_CropThreshCDCut.asc.gz",sep="")) & 
  #    !file.exists(paste(IndSppOutDir,"/",sp,"_","2090_Median_CropThreshCDCut.asc",sep=""))){
  #   file.remove(paste(IndSppOutDir,"/",sp,"_","2090_Median_CropThreshCDCut.asc",sep=""))
  # }
  # if(file.exists(paste(IndSppOutDir,"/",sp,"_","Current_CropThreshCDCut.asc.gz",sep="")) & 
  #    !file.exists(paste(IndSppOutDir,"/",sp,"_","Current_CropThreshCDCut.asc",sep=""))){
  #   file.remove(paste(IndSppOutDir,"/",sp,"_","Current_CropThreshCDCut.asc",sep=""))
  # }
  # 
  
  
  }
  #######################################################################
  RiskMed<-c(DatSumMod$Risk_Current[n]/DatSumMod$Risk_Current[n],
            DatSumMod$Risk_2050_Median[n]/DatSumMod$Risk_Current[n],
            DatSumMod$Risk_2090_Median[n]/DatSumMod$Risk_Current[n])
  #######################################################################
  ExpMed<-c(DatSumMod$Exposure_Current[n]/DatSumMod$Exposure_Current[n],
            DatSumMod$Exposure_2050_Median[n]/DatSumMod$Exposure_Current[n],
            DatSumMod$Exposure_2090_Median[n]/DatSumMod$Exposure_Current[n])
  ExpQ10<-c(DatSumMod$Exposure_Current[n]/DatSumMod$Exposure_Current[n],
            DatSumMod$Exposure_2050_Q10[n]/DatSumMod$Exposure_Current[n],
            DatSumMod$Exposure_2090_Q10[n]/DatSumMod$Exposure_Current[n])
  ExpQ90<-c(DatSumMod$Exposure_Current[n]/DatSumMod$Exposure_Current[n],
            DatSumMod$Exposure_2050_Q90[n]/DatSumMod$Exposure_Current[n],
            DatSumMod$Exposure_2090_Q90[n]/DatSumMod$Exposure_Current[n])
  #######################################################################
  ImpMed<-c(DatSumMod$Impact_Current[n]/DatSumMod$Impact_Current[n],
            DatSumMod$Impact_2050_Median[n]/DatSumMod$Impact_Current[n],
            DatSumMod$Impact_2090_Median[n]/DatSumMod$Impact_Current[n])

  
  
  #############################################
  
  ###########################################
  
  
  
  # png(filename=paste(PlotDir,"/RangeSize_CCTrends.png",sep=""), units="in", width=10, height=18, res=500) #png smaller file size, A3 for better resolution
  # par(mfrow = c(2,1), mar = c(1,1,3,1), oma=c(1,1,5,1), mgp=c(1.2,0.5,0))
  # #mar=c(bottom, left, top, right); oma= outer margins; mpg = axis label locations c(3, 1, 0)
  # 
  # 
  # 
  # #1. make species map
  # create vectors for species 2050 & 5090
  # map current expert distribution
  #plot 2050 & 2090 vectors to show shift

  
  #2. make species graph
  png(filename=paste(PlotDir,"/",sp,"_CCTrends.png",sep=""), units="in", width=14, height=7, res=500) #png smaller file size, A3 for better resolution
  par(mfrow = c(1,2), mar = c(3,3,0,0.5), oma=c(0.2,0.2,3,0.2), mgp=c(1.7,0.5,0))
  # #mar=c(bottom, left, top, right); oma= outer margins; mpg = axis label locations c(3, 1, 0)
  
  plot(myX,RangeMed, ylim=c(0,max(ExpQ90,SuitQ90,2,na.rm=TRUE)),xlim=c(1980,2090),col="white",bg="white", 
       bty="l", axes=TRUE, pch=21, lwd= 2, ylab="Fraction of current", xlab="year",lty=1, main=NULL)
  #plot zero line for reference
  lines (c(1980,2050,2090),c(1,1,1), col="grey", lwd= 2, lty=1)
  
  #2.a. plot % change in range size (black, solid thin)
  polygon(c(myX,rev(myX)),c(RangeQ10,rev(RangeQ90)), border = NA,
          col=rgb(col2rgb("black")[1],col2rgb("black")[2],col2rgb("black")[3],max=255,alpha=75))
  #points(myX,RangeMed, col = "black", bg="black",  pch=21, lwd= 1)
  lines (myX,RangeMed, col="black", lwd= 2, lty=1)
  lines (myX,RangeQ10, col="black", lwd= 1, lty=2)
  lines (myX,RangeQ90, col="black", lwd= 1, lty=2)
  
  #2.b. plot % of range lost (1985=0) (blue, solid thin)
  polygon(c(myX,rev(myX)),c(RLossQ10,rev(RLossQ90)), border = NA,
          col=rgb(col2rgb("royalblue")[1],col2rgb("royalblue")[2],col2rgb("royalblue")[3],max=255,alpha=75))
  #points(myX,RLossMed, col = "black", bg="royalblue",  pch=21, lwd= 1)
  lines (myX,RLossMed, col="royalblue", lwd= 2, lty=1)
  lines (myX,RLossQ10, col="royalblue", lwd= 1, lty=2)
  lines (myX,RLossQ90, col="royalblue", lwd= 1, lty=2)
  
  #2.c. plot % of range gained (1985=0) (red, solid thin)
  polygon(c(myX,rev(myX)),c(RGainQ10,rev(RGainQ90)),border = NA, 
          col=rgb(col2rgb("orangered")[1],col2rgb("orangered")[2],col2rgb("orangered")[3],max=255,alpha=75))
  #points(myX,RGainMed, col = "black", bg="orangered",  pch=21, lwd= 1)
  lines (myX,RGainMed, col="orangered", lwd= 2, lty=1)
  lines (myX,RGainQ10, col="orangered", lwd= 1, lty=2)
  lines (myX,RGainQ90, col="orangered", lwd= 1, lty=2)
  
  
  plot(myX,RangeMed, ylim=c(0,max(ExpQ90,SuitQ90,2,na.rm=TRUE)),xlim=c(1980,2090),col="white",bg="white", 
       bty="l", axes=TRUE, pch=21, ylab="Fraction of current", xlab="year", lwd= 2, lty=1, main=NULL)
  #plot zero line for reference
  lines (c(1970,2050,2100),c(1,1,1), col="grey", lwd= 1, lty=1)
  
  #2.d. plot % change in suitability/abundance (dashed)
  polygon(c(myX,rev(myX)),c(SuitQ10,rev(SuitQ90)), border = NA,
          col=rgb(col2rgb("darkgreen")[1],col2rgb("darkgreen")[2],col2rgb("darkgreen")[3],max=255,alpha=75))
  #points(myX,SuitMed, col = "black", bg="black",  pch=21, lwd= 1)
  lines (myX,SuitMed, col="darkgreen", lwd= 2, lty=1)
  lines (myX,SuitQ10, col="darkgreen", lwd= 1, lty=2)
  lines (myX,SuitQ90, col="darkgreen", lwd= 1, lty=2)
  #2.e. plot % of abundance lost (1985=0)?
  #2.f. plot % of abundance gained (1985=0)?
  
  #2.g. plot % change in human exposure (dotted)
  polygon(c(myX,rev(myX)),c(ExpQ10,rev(ExpQ90)), border = NA,
          col=rgb(col2rgb("orange")[1],col2rgb("orange")[2],col2rgb("orange")[3],max=255,alpha=75))
  #points(myX,ExpMed, col = "black", bg="black",  pch=21, lwd= 1)
  lines (myX,ExpMed, col="orange", lwd= 2, lty=1)
  lines (myX,ExpQ10, col="orange", lwd= 1, lty=2)
  lines (myX,ExpQ90, col="orange", lwd= 1, lty=2)
  #2.h. plot % new exposure (1985=0)?
  #2.i. plot % lost exposure (1985=0)?
  
  # #2.j. plot % change in impact (solid thick)
  # points(myX,ImpMed, col = "black", bg="black",  pch=23, lwd= 1)
  # lines (myX,ImpMed, col="black", lwd= 2, lty=1)
  # #2.k. plot % new impact (1985=0)
  # #2.l. plot % lost impact (1985=0)
  
  mtext(paste(gsub("_"," ",sp),sep=""), font= 4, side = 3, line = 1, outer = TRUE, cex= 1.5)
  
  dev.off()
  

    
  
  
  
  }

#####################
# 1. Range Size #####
#####################

{
png(filename=paste(PlotDir,"/RangeSize_CCTrends.png",sep=""), units="in", width=10, height=18, res=500) #png smaller file size, A3 for better resolution
par(mfrow = c(3,2), mar = c(1,1,3,1), oma=c(1,1,5,1), mgp=c(1.2,0.5,0))
#mar=c(bottom, left, top, right); oma= outer margins; mpg = axis label locations c(3, 1, 0)

myX<-c(1985,2050,2090)
DatSumMod$RangesizeSQKM_Current[DatSumMod$RangesizeSQKM_Current==0]<-0.01 #avoid dividing by 0 for % change graphs

# a. Median Range Size:
# a.i. raw

#plot empty base plot
myYMed<- sqrt((   c(DatSumMod$RangesizeSQKM_Current[length(DatSumMod[,1])],
                DatSumMod$RangesizeSQKM_2050_Median[length(DatSumMod[,1])],
                DatSumMod$RangesizeSQKM_2090_Median[length(DatSumMod[,1])])   )+1)
ymax<- sqrt(( max(max(DatSumMod$RangesizeSQKM_Current,na.rm=TRUE),
           max(DatSumMod$RangesizeSQKM_2050_Median,na.rm=TRUE),
           max(DatSumMod$RangesizeSQKM_2090_Median,na.rm=TRUE),na.rm=TRUE)     )+1)
plot(myX,myYMed, ylim=c(0,ymax),xlim=c(1970,2100),col="white",bg="white", type="b", 
     bty="n", axes=TRUE, pch=21, lwd= 2, lty=1, main="Median Range Size [sqrt(sqkm)]")

#plot each species' data
for (sp in spp){
  n<-which(DatSumMod$Accepted_Species==paste(gsub("_"," ",sp)))
  if(DatSumMod$RangesizeSQKM_Current[n]>0 & !is.na(DatSumMod$RangesizeSQKM_Current[n])){
    
    myYMed<-sqrt((  c(DatSumMod$RangesizeSQKM_Current[n],
                    DatSumMod$RangesizeSQKM_2050_Median[n],
                    DatSumMod$RangesizeSQKM_2090_Median[n])     )+1)
    
    points(myX,myYMed, col = "black", bg="darkgrey",  pch=21, lwd= 1)
    lines (myX,myYMed, col="darkgrey", lwd= 1, lty=1)
    
    if(!is.na(DatSumMod$RangesizeSQKM_2050_Median[n])){
      if(DatSumMod$RangesizeSQKM_2050_Median[n]>DatSumMod$RangesizeSQKM_Current[n]){
      lines (myX,myYMed, col="pink", lwd= 1, lty=1)
      print(paste(sp, " increase by 2050"))
    }}
    
    if(!is.na(DatSumMod$RangesizeSQKM_2090_Median[n])){
      if(DatSumMod$RangesizeSQKM_2090_Median[n]>DatSumMod$RangesizeSQKM_Current[n]){
      lines (myX,myYMed, col="red", lwd= 1, lty=1)
      print(paste(sp, " increase by 2090"))
      }}
    
}}

myYMed<-sqrt((   c(DatSumMod$RangesizeSQKM_Current[length(DatSumMod[,1])],
                 DatSumMod$RangesizeSQKM_2050_Median[length(DatSumMod[,1])],
                 DatSumMod$RangesizeSQKM_2090_Median[length(DatSumMod[,1])])    )+1)

points(myX,myYMed, col="black",bg="black",     pch=21, lwd= 2)
lines (myX,myYMed, col="black",                lwd= 2, lty=1)


# a.ii. percent of current

#plot empty base plot
myYMed<-c(DatSumMod$RangesizeSQKM_Current[length(DatSumMod[,1])]/DatSumMod$RangesizeSQKM_Current[length(DatSumMod[,1])],
          DatSumMod$RangesizeSQKM_2050_Median[length(DatSumMod[,1])]/DatSumMod$RangesizeSQKM_Current[length(DatSumMod[,1])],
          DatSumMod$RangesizeSQKM_2090_Median[length(DatSumMod[,1])]/DatSumMod$RangesizeSQKM_Current[length(DatSumMod[,1])])
ymax<-max(max(DatSumMod$RangesizeSQKM_Current/DatSumMod$RangesizeSQKM_Current,na.rm=TRUE),
          max(DatSumMod$RangesizeSQKM_2050_Median/DatSumMod$RangesizeSQKM_Current,na.rm=TRUE),
          max(DatSumMod$RangesizeSQKM_2090_Median/DatSumMod$RangesizeSQKM_Current,na.rm=TRUE),na.rm=TRUE)
plot(myX,myYMed, ylim=c(0,ymax),xlim=c(1970,2100),col="white",bg="white", type="b", 
     bty="n", axes=TRUE, pch=21, lwd= 2, lty=1, main="Median Range Size [% of Current]")

#plot each species' data
for (sp in spp){
  n<-which(DatSumMod$Accepted_Species==paste(gsub("_"," ",sp)))
  if(DatSumMod$RangesizeSQKM_Current[n]>0 & !is.na(DatSumMod$RangesizeSQKM_Current[n])){
    
    myYMed<-c(DatSumMod$RangesizeSQKM_Current[n]/DatSumMod$RangesizeSQKM_Current[n],
              DatSumMod$RangesizeSQKM_2050_Median[n]/DatSumMod$RangesizeSQKM_Current[n],
              DatSumMod$RangesizeSQKM_2090_Median[n]/DatSumMod$RangesizeSQKM_Current[n])
    
    points(myX,myYMed, col = "black", bg="darkgrey",  pch=21, lwd= 1)
    lines (myX,myYMed, col="darkgrey", lwd= 1, lty=1)
    
    if(!is.na(DatSumMod$RangesizeSQKM_2050_Median[n])){
      if(DatSumMod$RangesizeSQKM_2050_Median[n]>DatSumMod$RangesizeSQKM_Current[n]){
      lines (myX,myYMed, col="pink", lwd= 1, lty=1)
      print(paste(sp, " increase by 2050"))
    }}
    
    if(!is.na(DatSumMod$RangesizeSQKM_2090_Median[n])){
      if(DatSumMod$RangesizeSQKM_2090_Median[n]>DatSumMod$RangesizeSQKM_Current[n]){
      lines (myX,myYMed, col="red", lwd= 1, lty=1)
      print(paste(sp, " increase by 2090"))
      }}
    
  }}

myYMed<-c(DatSumMod$RangesizeSQKM_Current[length(DatSumMod[,1])]/DatSumMod$RangesizeSQKM_Current[length(DatSumMod[,1])],
          DatSumMod$RangesizeSQKM_2050_Median[length(DatSumMod[,1])]/DatSumMod$RangesizeSQKM_Current[length(DatSumMod[,1])],
          DatSumMod$RangesizeSQKM_2090_Median[length(DatSumMod[,1])]/DatSumMod$RangesizeSQKM_Current[length(DatSumMod[,1])])

points(myX,myYMed, col="black",bg="black",     pch=21, lwd= 2)
lines (myX,myYMed, col="black",                lwd= 2, lty=1)


# b. Q10 Range Size:
# b.i. raw

#plot empty base plot
myYMed<-sqrt((c(DatSumMod$RangesizeSQKM_Current[length(DatSumMod[,1])],
          DatSumMod$RangesizeSQKM_2050_Q10[length(DatSumMod[,1])],
          DatSumMod$RangesizeSQKM_2090_Q10[length(DatSumMod[,1])]) )+1)
ymax<-sqrt((max(max(DatSumMod$RangesizeSQKM_Current,na.rm=TRUE),
                max(DatSumMod$RangesizeSQKM_2050_Q10,na.rm=TRUE),
                max(DatSumMod$RangesizeSQKM_2090_Q10,na.rm=TRUE),na.rm=TRUE))+1)
plot(myX,myYMed, ylim=c(0,ymax),xlim=c(1970,2100),col="white",bg="white", type="b", 
     bty="n", axes=TRUE, pch=21, lwd= 2, lty=1, main="Q10 Range Size [sqrt(sqkm)]")

#plot each species' data
for (sp in spp){
  n<-which(DatSumMod$Accepted_Species==paste(gsub("_"," ",sp)))
  if(DatSumMod$RangesizeSQKM_Current[n]>0 & !is.na(DatSumMod$RangesizeSQKM_Current[n])){
    
    myYMed<-sqrt((c(DatSumMod$RangesizeSQKM_Current[n],
              DatSumMod$RangesizeSQKM_2050_Q10[n],
              DatSumMod$RangesizeSQKM_2090_Q10[n]))+1)
    
    points(myX,myYMed, col = "black", bg="darkgrey",  pch=21, lwd= 1)
    lines (myX,myYMed, col="darkgrey", lwd= 1, lty=1)
    
    if(!is.na(DatSumMod$RangesizeSQKM_2050_Q10[n])){
    if(DatSumMod$RangesizeSQKM_2050_Q10[n]>DatSumMod$RangesizeSQKM_Current[n]){
      lines (myX,myYMed, col="pink", lwd= 1, lty=1)
      print(paste(sp, " increase by 2050"))
    }}
    
    if(!is.na(DatSumMod$RangesizeSQKM_2090_Q10[n])){
      if(DatSumMod$RangesizeSQKM_2090_Q10[n]>DatSumMod$RangesizeSQKM_Current[n]){
      lines (myX,myYMed, col="red", lwd= 1, lty=1)
      print(paste(sp, " increase by 2090"))
      }}
    
  }}

myYMed<-sqrt((c(DatSumMod$RangesizeSQKM_Current[length(DatSumMod[,1])],
          DatSumMod$RangesizeSQKM_2050_Q10[length(DatSumMod[,1])],
          DatSumMod$RangesizeSQKM_2090_Q10[length(DatSumMod[,1])]))+1)

points(myX,myYMed, col="black",bg="black",     pch=21, lwd= 2)
lines (myX,myYMed, col="black",                lwd= 2, lty=1)


# b.ii. percent of current

#plot empty base plot
myYMed<-c(DatSumMod$RangesizeSQKM_Current[length(DatSumMod[,1])]/DatSumMod$RangesizeSQKM_Current[length(DatSumMod[,1])],
          DatSumMod$RangesizeSQKM_2050_Q10[length(DatSumMod[,1])]/DatSumMod$RangesizeSQKM_Current[length(DatSumMod[,1])],
          DatSumMod$RangesizeSQKM_2090_Q10[length(DatSumMod[,1])]/DatSumMod$RangesizeSQKM_Current[length(DatSumMod[,1])])
ymax<-max(max(DatSumMod$RangesizeSQKM_Current/DatSumMod$RangesizeSQKM_Current,na.rm=TRUE),
          max(DatSumMod$RangesizeSQKM_2050_Q10/DatSumMod$RangesizeSQKM_Current,na.rm=TRUE),
          max(DatSumMod$RangesizeSQKM_2090_Q10/DatSumMod$RangesizeSQKM_Current,na.rm=TRUE),na.rm=TRUE)
plot(myX,myYMed, ylim=c(0,ymax),xlim=c(1970,2100),col="white",bg="white", type="b", 
     bty="n", axes=TRUE, pch=21, lwd= 2, lty=1, main="Q10 Range Size [% of Current]")

#plot each species' data
for (sp in spp){
  n<-which(DatSumMod$Accepted_Species==paste(gsub("_"," ",sp)))
  if(DatSumMod$RangesizeSQKM_Current[n]>0 & !is.na(DatSumMod$RangesizeSQKM_Current[n])){
    
    myYMed<-c(DatSumMod$RangesizeSQKM_Current[n]/DatSumMod$RangesizeSQKM_Current[n],
              DatSumMod$RangesizeSQKM_2050_Q10[n]/DatSumMod$RangesizeSQKM_Current[n],
              DatSumMod$RangesizeSQKM_2090_Q10[n]/DatSumMod$RangesizeSQKM_Current[n])
    
    points(myX,myYMed, col = "black", bg="darkgrey",  pch=21, lwd= 1)
    lines (myX,myYMed, col="darkgrey", lwd= 1, lty=1)
    
    if(!is.na(DatSumMod$RangesizeSQKM_2050_Q10[n])){
      if(DatSumMod$RangesizeSQKM_2050_Q10[n]>DatSumMod$RangesizeSQKM_Current[n]){
      lines (myX,myYMed, col="pink", lwd= 1, lty=1)
      print(paste(sp, " increase by 2050"))
    }}
    
    if(!is.na(DatSumMod$RangesizeSQKM_2090_Q10[n])){
      if(DatSumMod$RangesizeSQKM_2090_Q10[n]>DatSumMod$RangesizeSQKM_Current[n]){
      lines (myX,myYMed, col="red", lwd= 1, lty=1)
      print(paste(sp, " increase by 2090"))
      
      }}
    
  }}

myYMed<-c(DatSumMod$RangesizeSQKM_Current[length(DatSumMod[,1])]/DatSumMod$RangesizeSQKM_Current[length(DatSumMod[,1])],
          DatSumMod$RangesizeSQKM_2050_Q10[length(DatSumMod[,1])]/DatSumMod$RangesizeSQKM_Current[length(DatSumMod[,1])],
          DatSumMod$RangesizeSQKM_2090_Q10[length(DatSumMod[,1])]/DatSumMod$RangesizeSQKM_Current[length(DatSumMod[,1])])

points(myX,myYMed, col="black",bg="black",     pch=21, lwd= 2)
lines (myX,myYMed, col="black",                lwd= 2, lty=1)

# c. Q90 Range Size:
# c.i. raw

#plot empty base plot
myYMed<-sqrt((c(DatSumMod$RangesizeSQKM_Current[length(DatSumMod[,1])],
          DatSumMod$RangesizeSQKM_2050_Q90[length(DatSumMod[,1])],
          DatSumMod$RangesizeSQKM_2090_Q90[length(DatSumMod[,1])]))+1)
ymax<-sqrt((max(max(DatSumMod$RangesizeSQKM_Current,na.rm=TRUE),
                max(DatSumMod$RangesizeSQKM_2050_Q90,na.rm=TRUE),
                max(DatSumMod$RangesizeSQKM_2090_Q90,na.rm=TRUE),na.rm=TRUE))+1)
plot(myX,myYMed, ylim=c(0,ymax),xlim=c(1970,2100),col="white",bg="white", type="b", 
     bty="n", axes=TRUE, pch=21, lwd= 2, lty=1, main="Q90 Range Size [sqrt(sqkm)]")

#plot each species' data
for (sp in spp){
  n<-which(DatSumMod$Accepted_Species==paste(gsub("_"," ",sp)))
  if(DatSumMod$RangesizeSQKM_Current[n]>0 & !is.na(DatSumMod$RangesizeSQKM_Current[n])){
    
    myYMed<-sqrt((c(DatSumMod$RangesizeSQKM_Current[n],
              DatSumMod$RangesizeSQKM_2050_Q90[n],
              DatSumMod$RangesizeSQKM_2090_Q90[n]))+1)
    
    points(myX,myYMed, col = "black", bg="darkgrey",  pch=21, lwd= 1)
    lines (myX,myYMed, col="darkgrey", lwd= 1, lty=1)
    
    if(!is.na(DatSumMod$RangesizeSQKM_2050_Q90[n])){
      if(DatSumMod$RangesizeSQKM_2050_Q90[n]>DatSumMod$RangesizeSQKM_Current[n]){
      lines (myX,myYMed, col="pink", lwd= 1, lty=1)
      print(paste(sp, " increase by 2050"))
    }}
    
    if(!is.na(DatSumMod$RangesizeSQKM_2090_Q90[n])){
      if(DatSumMod$RangesizeSQKM_2090_Q90[n]>DatSumMod$RangesizeSQKM_Current[n]){
      lines (myX,myYMed, col="red", lwd= 1, lty=1)
      print(paste(sp, " increase by 2090"))
      }}
    
  }}

myYMed<-sqrt((c(DatSumMod$RangesizeSQKM_Current[length(DatSumMod[,1])],
          DatSumMod$RangesizeSQKM_2050_Q90[length(DatSumMod[,1])],
          DatSumMod$RangesizeSQKM_2090_Q90[length(DatSumMod[,1])]))+1)

points(myX,myYMed, col="black",bg="black",     pch=21, lwd= 2)
lines (myX,myYMed, col="black",                lwd= 2, lty=1)


# c.ii. percent of current

#plot empty base plot
myYMed<-c(DatSumMod$RangesizeSQKM_Current[length(DatSumMod[,1])]/DatSumMod$RangesizeSQKM_Current[length(DatSumMod[,1])],
          DatSumMod$RangesizeSQKM_2050_Q90[length(DatSumMod[,1])]/DatSumMod$RangesizeSQKM_Current[length(DatSumMod[,1])],
          DatSumMod$RangesizeSQKM_2090_Q90[length(DatSumMod[,1])]/DatSumMod$RangesizeSQKM_Current[length(DatSumMod[,1])])
ymax<-max(max(DatSumMod$RangesizeSQKM_Current/DatSumMod$RangesizeSQKM_Current,na.rm=TRUE),
          max(DatSumMod$RangesizeSQKM_2050_Q90/DatSumMod$RangesizeSQKM_Current,na.rm=TRUE),
          max(DatSumMod$RangesizeSQKM_2090_Q90/DatSumMod$RangesizeSQKM_Current,na.rm=TRUE),na.rm=TRUE)
plot(myX,myYMed, ylim=c(0,ymax),xlim=c(1970,2100),col="white",bg="white", type="b", 
     bty="n", axes=TRUE, pch=21, lwd= 2, lty=1, main="Q90 Range Size [% of Current]")

#plot each species' data
for (sp in spp){
  n<-which(DatSumMod$Accepted_Species==paste(gsub("_"," ",sp)))
  if(DatSumMod$RangesizeSQKM_Current[n]>0 & !is.na(DatSumMod$RangesizeSQKM_Current[n])){
    
    myYMed<-c(DatSumMod$RangesizeSQKM_Current[n]/DatSumMod$RangesizeSQKM_Current[n],
              DatSumMod$RangesizeSQKM_2050_Q90[n]/DatSumMod$RangesizeSQKM_Current[n],
              DatSumMod$RangesizeSQKM_2090_Q90[n]/DatSumMod$RangesizeSQKM_Current[n])
    
    points(myX,myYMed, col = "black", bg="darkgrey",  pch=21, lwd= 1)
    lines (myX,myYMed, col="darkgrey", lwd= 1, lty=1)
    
    if(!is.na(DatSumMod$RangesizeSQKM_2050_Q90[n])){
      if(DatSumMod$RangesizeSQKM_2050_Q90[n]>DatSumMod$RangesizeSQKM_Current[n]){
      lines (myX,myYMed, col="pink", lwd= 1, lty=1)
      print(paste(sp, " increase by 2050"))
    }}
    
    if(!is.na(DatSumMod$RangesizeSQKM_2090_Q90[n])){
      if(DatSumMod$RangesizeSQKM_2090_Q90[n]>DatSumMod$RangesizeSQKM_Current[n]){
      lines (myX,myYMed, col="red", lwd= 1, lty=1)
      print(paste(sp, " increase by 2090"))
      }}
    
  }}

myYMed<-c(DatSumMod$RangesizeSQKM_Current[length(DatSumMod[,1])]/DatSumMod$RangesizeSQKM_Current[length(DatSumMod[,1])],
          DatSumMod$RangesizeSQKM_2050_Q90[length(DatSumMod[,1])]/DatSumMod$RangesizeSQKM_Current[length(DatSumMod[,1])],
          DatSumMod$RangesizeSQKM_2090_Q90[length(DatSumMod[,1])]/DatSumMod$RangesizeSQKM_Current[length(DatSumMod[,1])])

points(myX,myYMed, col="black",bg="black",     pch=21, lwd= 2)
lines (myX,myYMed, col="black",                lwd= 2, lty=1)

dev.off()

}



#####################
# 2. Suitability ####
#####################


{
  png(filename=paste(PlotDir,"/Suitability_CCTrends.png",sep=""), units="in", width=10, height=18, res=500) #png smaller file size, A3 for better resolution
  par(mfrow = c(3,2), mar = c(1,1,3,1), oma=c(1,1,5,1), mgp=c(1.2,0.5,0))
  #mar=c(bottom, left, top, right); oma= outer margins; mpg = axis label locations c(3, 1, 0)
  
  myX<-c(1985,2050,2090)
  DatSumMod$SuitabilitySum_Current[DatSumMod$SuitabilitySum_Current==0]<-0.01 #avoid dividing by 0 for % change graphs
  
  # a. Median Suitability:
  # a.i. raw
  
  #plot empty base plot
  myYMed<- sqrt((   c(DatSumMod$SuitabilitySum_Current[length(DatSumMod[,1])],
                      DatSumMod$SuitabilitySum_2050_Median[length(DatSumMod[,1])],
                      DatSumMod$SuitabilitySum_2090_Median[length(DatSumMod[,1])])   )+1)
  ymax<- sqrt(( max(max(DatSumMod$SuitabilitySum_Current,na.rm=TRUE),
                    max(DatSumMod$SuitabilitySum_2050_Median,na.rm=TRUE),
                    max(DatSumMod$SuitabilitySum_2090_Median,na.rm=TRUE),na.rm=TRUE)     )+1)
  plot(myX,myYMed, ylim=c(0,ymax),xlim=c(1970,2100),col="white",bg="white", type="b", 
       bty="n", axes=TRUE, pch=21, lwd= 2, lty=1, main="Median Suitability [sqrt(cell sum)]")
  
  #plot each species' data
  for (sp in spp){
    n<-which(DatSumMod$Accepted_Species==paste(gsub("_"," ",sp)))
    if(DatSumMod$SuitabilitySum_Current[n]>0 & !is.na(DatSumMod$SuitabilitySum_Current[n])){
      
      myYMed<-sqrt((  c(DatSumMod$SuitabilitySum_Current[n],
                        DatSumMod$SuitabilitySum_2050_Median[n],
                        DatSumMod$SuitabilitySum_2090_Median[n])     )+1)
      
      points(myX,myYMed, col = "black", bg="darkgrey",  pch=21, lwd= 1)
      lines (myX,myYMed, col="darkgrey", lwd= 1, lty=1)
      
      if(!is.na(DatSumMod$SuitabilitySum_2050_Median[n])){
        if(DatSumMod$SuitabilitySum_2050_Median[n]>DatSumMod$SuitabilitySum_Current[n]){
          lines (myX,myYMed, col="pink", lwd= 1, lty=1)
          print(paste(sp, " increase by 2050"))
        }}
      
      if(!is.na(DatSumMod$SuitabilitySum_2090_Median[n])){
        if(DatSumMod$SuitabilitySum_2090_Median[n]>DatSumMod$SuitabilitySum_Current[n]){
          lines (myX,myYMed, col="red", lwd= 1, lty=1)
          print(paste(sp, " increase by 2090"))
        }}
      
    }}
  
  myYMed<-sqrt((   c(DatSumMod$SuitabilitySum_Current[length(DatSumMod[,1])],
                     DatSumMod$SuitabilitySum_2050_Median[length(DatSumMod[,1])],
                     DatSumMod$SuitabilitySum_2090_Median[length(DatSumMod[,1])])    )+1)
  
  points(myX,myYMed, col="black",bg="black",     pch=21, lwd= 2)
  lines (myX,myYMed, col="black",                lwd= 2, lty=1)
  
  
  # a.ii. percent of current
  
  #plot empty base plot
  myYMed<-c(DatSumMod$SuitabilitySum_Current[length(DatSumMod[,1])]/DatSumMod$SuitabilitySum_Current[length(DatSumMod[,1])],
            DatSumMod$SuitabilitySum_2050_Median[length(DatSumMod[,1])]/DatSumMod$SuitabilitySum_Current[length(DatSumMod[,1])],
            DatSumMod$SuitabilitySum_2090_Median[length(DatSumMod[,1])]/DatSumMod$SuitabilitySum_Current[length(DatSumMod[,1])])
  ymax<-max(max(DatSumMod$SuitabilitySum_Current/DatSumMod$SuitabilitySum_Current,na.rm=TRUE),
            max(DatSumMod$SuitabilitySum_2050_Median/DatSumMod$SuitabilitySum_Current,na.rm=TRUE),
            max(DatSumMod$SuitabilitySum_2090_Median/DatSumMod$SuitabilitySum_Current,na.rm=TRUE),na.rm=TRUE)
  plot(myX,myYMed, ylim=c(0,ymax),xlim=c(1970,2100),col="white",bg="white", type="b", 
       bty="n", axes=TRUE, pch=21, lwd= 2, lty=1, main="Median Suitability [% of Current]")
  
  #plot each species' data
  for (sp in spp){
    n<-which(DatSumMod$Accepted_Species==paste(gsub("_"," ",sp)))
    if(DatSumMod$SuitabilitySum_Current[n]>0 & !is.na(DatSumMod$SuitabilitySum_Current[n])){
      
      myYMed<-c(DatSumMod$SuitabilitySum_Current[n]/DatSumMod$SuitabilitySum_Current[n],
                DatSumMod$SuitabilitySum_2050_Median[n]/DatSumMod$SuitabilitySum_Current[n],
                DatSumMod$SuitabilitySum_2090_Median[n]/DatSumMod$SuitabilitySum_Current[n])
      
      points(myX,myYMed, col = "black", bg="darkgrey",  pch=21, lwd= 1)
      lines (myX,myYMed, col="darkgrey", lwd= 1, lty=1)
      
      if(!is.na(DatSumMod$SuitabilitySum_2050_Median[n])){
        if(DatSumMod$SuitabilitySum_2050_Median[n]>DatSumMod$SuitabilitySum_Current[n]){
          lines (myX,myYMed, col="pink", lwd= 1, lty=1)
          print(paste(sp, " increase by 2050"))
        }}
      
      if(!is.na(DatSumMod$SuitabilitySum_2090_Median[n])){
        if(DatSumMod$SuitabilitySum_2090_Median[n]>DatSumMod$SuitabilitySum_Current[n]){
          lines (myX,myYMed, col="red", lwd= 1, lty=1)
          print(paste(sp, " increase by 2090"))
        }}
      
    }}
  
  myYMed<-c(DatSumMod$SuitabilitySum_Current[length(DatSumMod[,1])]/DatSumMod$SuitabilitySum_Current[length(DatSumMod[,1])],
            DatSumMod$SuitabilitySum_2050_Median[length(DatSumMod[,1])]/DatSumMod$SuitabilitySum_Current[length(DatSumMod[,1])],
            DatSumMod$SuitabilitySum_2090_Median[length(DatSumMod[,1])]/DatSumMod$SuitabilitySum_Current[length(DatSumMod[,1])])
  
  points(myX,myYMed, col="black",bg="black",     pch=21, lwd= 2)
  lines (myX,myYMed, col="black",                lwd= 2, lty=1)
  
  
  # b. Q10 Suitability:
  # b.i. raw
  
  #plot empty base plot
  myYMed<-sqrt((c(DatSumMod$SuitabilitySum_Current[length(DatSumMod[,1])],
                  DatSumMod$SuitabilitySum_2050_Q10[length(DatSumMod[,1])],
                  DatSumMod$SuitabilitySum_2090_Q10[length(DatSumMod[,1])]) )+1)
  ymax<-sqrt((max(max(DatSumMod$SuitabilitySum_Current,na.rm=TRUE),
                  max(DatSumMod$SuitabilitySum_2050_Q10,na.rm=TRUE),
                  max(DatSumMod$SuitabilitySum_2090_Q10,na.rm=TRUE),na.rm=TRUE))+1)
  plot(myX,myYMed, ylim=c(0,ymax),xlim=c(1970,2100),col="white",bg="white", type="b", 
       bty="n", axes=TRUE, pch=21, lwd= 2, lty=1, main="Q10 Suitability [sqrt(cell sum)]")
  
  #plot each species' data
  for (sp in spp){
    n<-which(DatSumMod$Accepted_Species==paste(gsub("_"," ",sp)))
    if(DatSumMod$SuitabilitySum_Current[n]>0 & !is.na(DatSumMod$SuitabilitySum_Current[n])){
      
      myYMed<-sqrt((c(DatSumMod$SuitabilitySum_Current[n],
                      DatSumMod$SuitabilitySum_2050_Q10[n],
                      DatSumMod$SuitabilitySum_2090_Q10[n]))+1)
      
      points(myX,myYMed, col = "black", bg="darkgrey",  pch=21, lwd= 1)
      lines (myX,myYMed, col="darkgrey", lwd= 1, lty=1)
      
      if(!is.na(DatSumMod$SuitabilitySum_2050_Q10[n])){
        if(DatSumMod$SuitabilitySum_2050_Q10[n]>DatSumMod$SuitabilitySum_Current[n]){
          lines (myX,myYMed, col="pink", lwd= 1, lty=1)
          print(paste(sp, " increase by 2050"))
        }}
      
      if(!is.na(DatSumMod$SuitabilitySum_2090_Q10[n])){
        if(DatSumMod$SuitabilitySum_2090_Q10[n]>DatSumMod$SuitabilitySum_Current[n]){
          lines (myX,myYMed, col="red", lwd= 1, lty=1)
          print(paste(sp, " increase by 2090"))
        }}
      
    }}
  
  myYMed<-sqrt((c(DatSumMod$SuitabilitySum_Current[length(DatSumMod[,1])],
                  DatSumMod$SuitabilitySum_2050_Q10[length(DatSumMod[,1])],
                  DatSumMod$SuitabilitySum_2090_Q10[length(DatSumMod[,1])]))+1)
  
  points(myX,myYMed, col="black",bg="black",     pch=21, lwd= 2)
  lines (myX,myYMed, col="black",                lwd= 2, lty=1)
  
  
  # b.ii. percent of current
  
  #plot empty base plot
  myYMed<-c(DatSumMod$SuitabilitySum_Current[length(DatSumMod[,1])]/DatSumMod$SuitabilitySum_Current[length(DatSumMod[,1])],
            DatSumMod$SuitabilitySum_2050_Q10[length(DatSumMod[,1])]/DatSumMod$SuitabilitySum_Current[length(DatSumMod[,1])],
            DatSumMod$SuitabilitySum_2090_Q10[length(DatSumMod[,1])]/DatSumMod$SuitabilitySum_Current[length(DatSumMod[,1])])
  ymax<-max(max(DatSumMod$SuitabilitySum_Current/DatSumMod$SuitabilitySum_Current,na.rm=TRUE),
            max(DatSumMod$SuitabilitySum_2050_Q10/DatSumMod$SuitabilitySum_Current,na.rm=TRUE),
            max(DatSumMod$SuitabilitySum_2090_Q10/DatSumMod$SuitabilitySum_Current,na.rm=TRUE),na.rm=TRUE)
  plot(myX,myYMed, ylim=c(0,ymax),xlim=c(1970,2100),col="white",bg="white", type="b", 
       bty="n", axes=TRUE, pch=21, lwd= 2, lty=1, main="Q10 Suitability [% of Current]")
  
  #plot each species' data
  for (sp in spp){
    n<-which(DatSumMod$Accepted_Species==paste(gsub("_"," ",sp)))
    if(DatSumMod$SuitabilitySum_Current[n]>0 & !is.na(DatSumMod$SuitabilitySum_Current[n])){
      
      myYMed<-c(DatSumMod$SuitabilitySum_Current[n]/DatSumMod$SuitabilitySum_Current[n],
                DatSumMod$SuitabilitySum_2050_Q10[n]/DatSumMod$SuitabilitySum_Current[n],
                DatSumMod$SuitabilitySum_2090_Q10[n]/DatSumMod$SuitabilitySum_Current[n])
      
      points(myX,myYMed, col = "black", bg="darkgrey",  pch=21, lwd= 1)
      lines (myX,myYMed, col="darkgrey", lwd= 1, lty=1)
      
      if(!is.na(DatSumMod$SuitabilitySum_2050_Q10[n])){
        if(DatSumMod$SuitabilitySum_2050_Q10[n]>DatSumMod$SuitabilitySum_Current[n]){
          lines (myX,myYMed, col="pink", lwd= 1, lty=1)
          print(paste(sp, " increase by 2050"))
        }}
      
      if(!is.na(DatSumMod$SuitabilitySum_2090_Q10[n])){
        if(DatSumMod$SuitabilitySum_2090_Q10[n]>DatSumMod$SuitabilitySum_Current[n]){
          lines (myX,myYMed, col="red", lwd= 1, lty=1)
          print(paste(sp, " increase by 2090"))
          
        }}
      
    }}
  
  myYMed<-c(DatSumMod$SuitabilitySum_Current[length(DatSumMod[,1])]/DatSumMod$SuitabilitySum_Current[length(DatSumMod[,1])],
            DatSumMod$SuitabilitySum_2050_Q10[length(DatSumMod[,1])]/DatSumMod$SuitabilitySum_Current[length(DatSumMod[,1])],
            DatSumMod$SuitabilitySum_2090_Q10[length(DatSumMod[,1])]/DatSumMod$SuitabilitySum_Current[length(DatSumMod[,1])])
  
  points(myX,myYMed, col="black",bg="black",     pch=21, lwd= 2)
  lines (myX,myYMed, col="black",                lwd= 2, lty=1)
  
  # c. Q90 Suitability:
  # c.i. raw
  
  #plot empty base plot
  myYMed<-sqrt((c(DatSumMod$SuitabilitySum_Current[length(DatSumMod[,1])],
                  DatSumMod$SuitabilitySum_2050_Q90[length(DatSumMod[,1])],
                  DatSumMod$SuitabilitySum_2090_Q90[length(DatSumMod[,1])]))+1)
  ymax<-sqrt((max(max(DatSumMod$SuitabilitySum_Current,na.rm=TRUE),
                  max(DatSumMod$SuitabilitySum_2050_Q90,na.rm=TRUE),
                  max(DatSumMod$SuitabilitySum_2090_Q90,na.rm=TRUE),na.rm=TRUE))+1)
  plot(myX,myYMed, ylim=c(0,ymax),xlim=c(1970,2100),col="white",bg="white", type="b", 
       bty="n", axes=TRUE, pch=21, lwd= 2, lty=1, main="Q90 SuitabilitySum [sqrt(cell sum)]")
  
  #plot each species' data
  for (sp in spp){
    n<-which(DatSumMod$Accepted_Species==paste(gsub("_"," ",sp)))
    if(DatSumMod$SuitabilitySum_Current[n]>0 & !is.na(DatSumMod$SuitabilitySum_Current[n])){
      
      myYMed<-sqrt((c(DatSumMod$SuitabilitySum_Current[n],
                      DatSumMod$SuitabilitySum_2050_Q90[n],
                      DatSumMod$SuitabilitySum_2090_Q90[n]))+1)
      
      points(myX,myYMed, col = "black", bg="darkgrey",  pch=21, lwd= 1)
      lines (myX,myYMed, col="darkgrey", lwd= 1, lty=1)
      
      if(!is.na(DatSumMod$SuitabilitySum_2050_Q90[n])){
        if(DatSumMod$SuitabilitySum_2050_Q90[n]>DatSumMod$SuitabilitySum_Current[n]){
          lines (myX,myYMed, col="pink", lwd= 1, lty=1)
          print(paste(sp, " increase by 2050"))
        }}
      
      if(!is.na(DatSumMod$SuitabilitySum_2090_Q90[n])){
        if(DatSumMod$SuitabilitySum_2090_Q90[n]>DatSumMod$SuitabilitySum_Current[n]){
          lines (myX,myYMed, col="red", lwd= 1, lty=1)
          print(paste(sp, " increase by 2090"))
        }}
      
    }}
  
  myYMed<-sqrt((c(DatSumMod$SuitabilitySum_Current[length(DatSumMod[,1])],
                  DatSumMod$SuitabilitySum_2050_Q90[length(DatSumMod[,1])],
                  DatSumMod$SuitabilitySum_2090_Q90[length(DatSumMod[,1])]))+1)
  
  points(myX,myYMed, col="black",bg="black",     pch=21, lwd= 2)
  lines (myX,myYMed, col="black",                lwd= 2, lty=1)
  
  
  # c.ii. percent of current
  
  #plot empty base plot
  myYMed<-c(DatSumMod$SuitabilitySum_Current[length(DatSumMod[,1])]/DatSumMod$SuitabilitySum_Current[length(DatSumMod[,1])],
            DatSumMod$SuitabilitySum_2050_Q90[length(DatSumMod[,1])]/DatSumMod$SuitabilitySum_Current[length(DatSumMod[,1])],
            DatSumMod$SuitabilitySum_2090_Q90[length(DatSumMod[,1])]/DatSumMod$SuitabilitySum_Current[length(DatSumMod[,1])])
  ymax<-max(max(DatSumMod$SuitabilitySum_Current/DatSumMod$SuitabilitySum_Current,na.rm=TRUE),
            max(DatSumMod$SuitabilitySum_2050_Q90/DatSumMod$SuitabilitySum_Current,na.rm=TRUE),
            max(DatSumMod$SuitabilitySum_2090_Q90/DatSumMod$SuitabilitySum_Current,na.rm=TRUE),na.rm=TRUE)
  plot(myX,myYMed, ylim=c(0,ymax),xlim=c(1970,2100),col="white",bg="white", type="b", 
       bty="n", axes=TRUE, pch=21, lwd= 2, lty=1, main="Q90 Suitability [% of Current]")
  
  #plot each species' data
  for (sp in spp){
    n<-which(DatSumMod$Accepted_Species==paste(gsub("_"," ",sp)))
    if(DatSumMod$SuitabilitySum_Current[n]>0 & !is.na(DatSumMod$SuitabilitySum_Current[n])){
      
      myYMed<-c(DatSumMod$SuitabilitySum_Current[n]/DatSumMod$SuitabilitySum_Current[n],
                DatSumMod$SuitabilitySum_2050_Q90[n]/DatSumMod$SuitabilitySum_Current[n],
                DatSumMod$SuitabilitySum_2090_Q90[n]/DatSumMod$SuitabilitySum_Current[n])
      
      points(myX,myYMed, col = "black", bg="darkgrey",  pch=21, lwd= 1)
      lines (myX,myYMed, col="darkgrey", lwd= 1, lty=1)
      
      if(!is.na(DatSumMod$SuitabilitySum_2050_Q90[n])){
        if(DatSumMod$SuitabilitySum_2050_Q90[n]>DatSumMod$SuitabilitySum_Current[n]){
          lines (myX,myYMed, col="pink", lwd= 1, lty=1)
          print(paste(sp, " increase by 2050"))
        }}
      
      if(!is.na(DatSumMod$SuitabilitySum_2090_Q90[n])){
        if(DatSumMod$SuitabilitySum_2090_Q90[n]>DatSumMod$SuitabilitySum_Current[n]){
          lines (myX,myYMed, col="red", lwd= 1, lty=1)
          print(paste(sp, " increase by 2090"))
        }}
      
    }}
  
  myYMed<-c(DatSumMod$SuitabilitySum_Current[length(DatSumMod[,1])]/DatSumMod$SuitabilitySum_Current[length(DatSumMod[,1])],
            DatSumMod$SuitabilitySum_2050_Q90[length(DatSumMod[,1])]/DatSumMod$SuitabilitySum_Current[length(DatSumMod[,1])],
            DatSumMod$SuitabilitySum_2090_Q90[length(DatSumMod[,1])]/DatSumMod$SuitabilitySum_Current[length(DatSumMod[,1])])
  
  points(myX,myYMed, col="black",bg="black",     pch=21, lwd= 2)
  lines (myX,myYMed, col="black",                lwd= 2, lty=1)
  
  dev.off()
  
}


#####################
# 3. Exposure #######
#####################


{
  png(filename=paste(PlotDir,"/Exposure_CCTrends.png",sep=""), units="in", width=10, height=18, res=500) #png smaller file size, A3 for better resolution
  par(mfrow = c(3,2), mar = c(1,1,3,1), oma=c(1,1,5,1), mgp=c(1.2,0.5,0))
  #mar=c(bottom, left, top, right); oma= outer margins; mpg = axis label locations c(3, 1, 0)
  
  myX<-c(1985,2050,2090)
  DatSumMod$Exposure_Current[DatSumMod$Exposure_Current==0]<-0.01 #avoid dividing by 0 for % change graphs
  
  # a. Median Exposure:
  # a.i. raw
  
  #plot empty base plot
  myYMed<- sqrt((   c(DatSumMod$Exposure_Current[length(DatSumMod[,1])],
                      DatSumMod$Exposure_2050_Median[length(DatSumMod[,1])],
                      DatSumMod$Exposure_2090_Median[length(DatSumMod[,1])])   )+1)
  ymax<- sqrt(( max(max(DatSumMod$Exposure_Current,na.rm=TRUE),
                    max(DatSumMod$Exposure_2050_Median,na.rm=TRUE),
                    max(DatSumMod$Exposure_2090_Median,na.rm=TRUE),na.rm=TRUE)     )+1)
  plot(myX,myYMed, ylim=c(0,ymax),xlim=c(1970,2100),col="white",bg="white", type="b", 
       bty="n", axes=TRUE, pch=21, lwd= 2, lty=1, main="Median Exposure [sqrt(cell sum)]")
  
  #plot each species' data
  for (sp in spp){
    n<-which(DatSumMod$Accepted_Species==paste(gsub("_"," ",sp)))
    if(DatSumMod$Exposure_Current[n]>0 & !is.na(DatSumMod$Exposure_Current[n])){
      
      myYMed<-sqrt((  c(DatSumMod$Exposure_Current[n],
                        DatSumMod$Exposure_2050_Median[n],
                        DatSumMod$Exposure_2090_Median[n])     )+1)
      
      points(myX,myYMed, col = "black", bg="darkgrey",  pch=21, lwd= 1)
      lines (myX,myYMed, col="darkgrey", lwd= 1, lty=1)
      
      if(!is.na(DatSumMod$Exposure_2050_Median[n])){
        if(DatSumMod$Exposure_2050_Median[n]>DatSumMod$Exposure_Current[n]){
          lines (myX,myYMed, col="pink", lwd= 1, lty=1)
          print(paste(sp, " increase by 2050"))
        }}
      
      if(!is.na(DatSumMod$Exposure_2090_Median[n])){
        if(DatSumMod$Exposure_2090_Median[n]>DatSumMod$Exposure_Current[n]){
          lines (myX,myYMed, col="red", lwd= 1, lty=1)
          print(paste(sp, " increase by 2090"))
        }}
      
    }}
  
  myYMed<-sqrt((   c(DatSumMod$Exposure_Current[length(DatSumMod[,1])],
                     DatSumMod$Exposure_2050_Median[length(DatSumMod[,1])],
                     DatSumMod$Exposure_2090_Median[length(DatSumMod[,1])])    )+1)
  
  points(myX,myYMed, col="black",bg="black",     pch=21, lwd= 2)
  lines (myX,myYMed, col="black",                lwd= 2, lty=1)
  
  
  # a.ii. percent of current
  
  #plot empty base plot
  myYMed<-c(DatSumMod$Exposure_Current[length(DatSumMod[,1])]/DatSumMod$Exposure_Current[length(DatSumMod[,1])],
            DatSumMod$Exposure_2050_Median[length(DatSumMod[,1])]/DatSumMod$Exposure_Current[length(DatSumMod[,1])],
            DatSumMod$Exposure_2090_Median[length(DatSumMod[,1])]/DatSumMod$Exposure_Current[length(DatSumMod[,1])])
  ymax<-max(max(DatSumMod$Exposure_Current/DatSumMod$Exposure_Current,na.rm=TRUE),
            max(DatSumMod$Exposure_2050_Median/DatSumMod$Exposure_Current,na.rm=TRUE),
            max(DatSumMod$Exposure_2090_Median/DatSumMod$Exposure_Current,na.rm=TRUE),na.rm=TRUE)
  plot(myX,myYMed, ylim=c(0,ymax),xlim=c(1970,2100),col="white",bg="white", type="b", 
       bty="n", axes=TRUE, pch=21, lwd= 2, lty=1, main="Median Exposure [% of Current]")
  
  #plot each species' data
  for (sp in spp){
    n<-which(DatSumMod$Accepted_Species==paste(gsub("_"," ",sp)))
    if(DatSumMod$Exposure_Current[n]>0 & !is.na(DatSumMod$Exposure_Current[n])){
      
      myYMed<-c(DatSumMod$Exposure_Current[n]/DatSumMod$Exposure_Current[n],
                DatSumMod$Exposure_2050_Median[n]/DatSumMod$Exposure_Current[n],
                DatSumMod$Exposure_2090_Median[n]/DatSumMod$Exposure_Current[n])
      
      points(myX,myYMed, col = "black", bg="darkgrey",  pch=21, lwd= 1)
      lines (myX,myYMed, col="darkgrey", lwd= 1, lty=1)
      
      if(!is.na(DatSumMod$Exposure_2050_Median[n])){
        if(DatSumMod$Exposure_2050_Median[n]>DatSumMod$Exposure_Current[n]){
          lines (myX,myYMed, col="pink", lwd= 1, lty=1)
          print(paste(sp, " increase by 2050"))
        }}
      
      if(!is.na(DatSumMod$Exposure_2090_Median[n])){
        if(DatSumMod$Exposure_2090_Median[n]>DatSumMod$Exposure_Current[n]){
          lines (myX,myYMed, col="red", lwd= 1, lty=1)
          print(paste(sp, " increase by 2090"))
        }}
      
    }}
  
  myYMed<-c(DatSumMod$Exposure_Current[length(DatSumMod[,1])]/DatSumMod$Exposure_Current[length(DatSumMod[,1])],
            DatSumMod$Exposure_2050_Median[length(DatSumMod[,1])]/DatSumMod$Exposure_Current[length(DatSumMod[,1])],
            DatSumMod$Exposure_2090_Median[length(DatSumMod[,1])]/DatSumMod$Exposure_Current[length(DatSumMod[,1])])
  
  points(myX,myYMed, col="black",bg="black",     pch=21, lwd= 2)
  lines (myX,myYMed, col="black",                lwd= 2, lty=1)
  
  
  # b. Q10 Exposure:
  # b.i. raw
  
  #plot empty base plot
  myYMed<-sqrt((c(DatSumMod$Exposure_Current[length(DatSumMod[,1])],
                  DatSumMod$Exposure_2050_Q10[length(DatSumMod[,1])],
                  DatSumMod$Exposure_2090_Q10[length(DatSumMod[,1])]) )+1)
  ymax<-sqrt((max(max(DatSumMod$Exposure_Current,na.rm=TRUE),
                  max(DatSumMod$Exposure_2050_Q10,na.rm=TRUE),
                  max(DatSumMod$Exposure_2090_Q10,na.rm=TRUE),na.rm=TRUE))+1)
  plot(myX,myYMed, ylim=c(0,ymax),xlim=c(1970,2100),col="white",bg="white", type="b", 
       bty="n", axes=TRUE, pch=21, lwd= 2, lty=1, main="Q10 Exposure [sqrt(cell sum)]")
  
  #plot each species' data
  for (sp in spp){
    n<-which(DatSumMod$Accepted_Species==paste(gsub("_"," ",sp)))
    if(DatSumMod$Exposure_Current[n]>0 & !is.na(DatSumMod$Exposure_Current[n])){
      
      myYMed<-sqrt((c(DatSumMod$Exposure_Current[n],
                      DatSumMod$Exposure_2050_Q10[n],
                      DatSumMod$Exposure_2090_Q10[n]))+1)
      
      points(myX,myYMed, col = "black", bg="darkgrey",  pch=21, lwd= 1)
      lines (myX,myYMed, col="darkgrey", lwd= 1, lty=1)
      
      if(!is.na(DatSumMod$Exposure_2050_Q10[n])){
        if(DatSumMod$Exposure_2050_Q10[n]>DatSumMod$Exposure_Current[n]){
          lines (myX,myYMed, col="pink", lwd= 1, lty=1)
          print(paste(sp, " increase by 2050"))
        }}
      
      if(!is.na(DatSumMod$Exposure_2090_Q10[n])){
        if(DatSumMod$Exposure_2090_Q10[n]>DatSumMod$Exposure_Current[n]){
          lines (myX,myYMed, col="red", lwd= 1, lty=1)
          print(paste(sp, " increase by 2090"))
        }}
      
    }}
  
  myYMed<-sqrt((c(DatSumMod$Exposure_Current[length(DatSumMod[,1])],
                  DatSumMod$Exposure_2050_Q10[length(DatSumMod[,1])],
                  DatSumMod$Exposure_2090_Q10[length(DatSumMod[,1])]))+1)
  
  points(myX,myYMed, col="black",bg="black",     pch=21, lwd= 2)
  lines (myX,myYMed, col="black",                lwd= 2, lty=1)
  
  
  # b.ii. percent of current
  
  #plot empty base plot
  myYMed<-c(DatSumMod$Exposure_Current[length(DatSumMod[,1])]/DatSumMod$Exposure_Current[length(DatSumMod[,1])],
            DatSumMod$Exposure_2050_Q10[length(DatSumMod[,1])]/DatSumMod$Exposure_Current[length(DatSumMod[,1])],
            DatSumMod$Exposure_2090_Q10[length(DatSumMod[,1])]/DatSumMod$Exposure_Current[length(DatSumMod[,1])])
  ymax<-max(max(DatSumMod$Exposure_Current/DatSumMod$Exposure_Current,na.rm=TRUE),
            max(DatSumMod$Exposure_2050_Q10/DatSumMod$Exposure_Current,na.rm=TRUE),
            max(DatSumMod$Exposure_2090_Q10/DatSumMod$Exposure_Current,na.rm=TRUE),na.rm=TRUE)
  plot(myX,myYMed, ylim=c(0,ymax),xlim=c(1970,2100),col="white",bg="white", type="b", 
       bty="n", axes=TRUE, pch=21, lwd= 2, lty=1, main="Q10 Exposure [% of Current]")
  
  #plot each species' data
  for (sp in spp){
    n<-which(DatSumMod$Accepted_Species==paste(gsub("_"," ",sp)))
    if(DatSumMod$Exposure_Current[n]>0 & !is.na(DatSumMod$Exposure_Current[n])){
      
      myYMed<-c(DatSumMod$Exposure_Current[n]/DatSumMod$Exposure_Current[n],
                DatSumMod$Exposure_2050_Q10[n]/DatSumMod$Exposure_Current[n],
                DatSumMod$Exposure_2090_Q10[n]/DatSumMod$Exposure_Current[n])
      
      points(myX,myYMed, col = "black", bg="darkgrey",  pch=21, lwd= 1)
      lines (myX,myYMed, col="darkgrey", lwd= 1, lty=1)
      
      if(!is.na(DatSumMod$Exposure_2050_Q10[n])){
        if(DatSumMod$Exposure_2050_Q10[n]>DatSumMod$Exposure_Current[n]){
          lines (myX,myYMed, col="pink", lwd= 1, lty=1)
          print(paste(sp, " increase by 2050"))
        }}
      
      if(!is.na(DatSumMod$Exposure_2090_Q10[n])){
        if(DatSumMod$Exposure_2090_Q10[n]>DatSumMod$Exposure_Current[n]){
          lines (myX,myYMed, col="red", lwd= 1, lty=1)
          print(paste(sp, " increase by 2090"))
          
        }}
      
    }}
  
  myYMed<-c(DatSumMod$Exposure_Current[length(DatSumMod[,1])]/DatSumMod$Exposure_Current[length(DatSumMod[,1])],
            DatSumMod$Exposure_2050_Q10[length(DatSumMod[,1])]/DatSumMod$Exposure_Current[length(DatSumMod[,1])],
            DatSumMod$Exposure_2090_Q10[length(DatSumMod[,1])]/DatSumMod$Exposure_Current[length(DatSumMod[,1])])
  
  points(myX,myYMed, col="black",bg="black",     pch=21, lwd= 2)
  lines (myX,myYMed, col="black",                lwd= 2, lty=1)
  
  # c. Q90 Exposure:
  # c.i. raw
  
  #plot empty base plot
  myYMed<-sqrt((c(DatSumMod$Exposure_Current[length(DatSumMod[,1])],
                  DatSumMod$Exposure_2050_Q90[length(DatSumMod[,1])],
                  DatSumMod$Exposure_2090_Q90[length(DatSumMod[,1])]))+1)
  ymax<-sqrt((max(max(DatSumMod$Exposure_Current,na.rm=TRUE),
                  max(DatSumMod$Exposure_2050_Q90,na.rm=TRUE),
                  max(DatSumMod$Exposure_2090_Q90,na.rm=TRUE),na.rm=TRUE))+1)
  plot(myX,myYMed, ylim=c(0,ymax),xlim=c(1970,2100),col="white",bg="white", type="b", 
       bty="n", axes=TRUE, pch=21, lwd= 2, lty=1, main="Q90 Exposure [sqrt(cell sum)]")
  
  #plot each species' data
  for (sp in spp){
    n<-which(DatSumMod$Accepted_Species==paste(gsub("_"," ",sp)))
    if(DatSumMod$Exposure_Current[n]>0 & !is.na(DatSumMod$Exposure_Current[n])){
      
      myYMed<-sqrt((c(DatSumMod$Exposure_Current[n],
                      DatSumMod$Exposure_2050_Q90[n],
                      DatSumMod$Exposure_2090_Q90[n]))+1)
      
      points(myX,myYMed, col = "black", bg="darkgrey",  pch=21, lwd= 1)
      lines (myX,myYMed, col="darkgrey", lwd= 1, lty=1)
      
      if(!is.na(DatSumMod$Exposure_2050_Q90[n])){
        if(DatSumMod$Exposure_2050_Q90[n]>DatSumMod$Exposure_Current[n]){
          lines (myX,myYMed, col="pink", lwd= 1, lty=1)
          print(paste(sp, " increase by 2050"))
        }}
      
      if(!is.na(DatSumMod$Exposure_2090_Q90[n])){
        if(DatSumMod$Exposure_2090_Q90[n]>DatSumMod$Exposure_Current[n]){
          lines (myX,myYMed, col="red", lwd= 1, lty=1)
          print(paste(sp, " increase by 2090"))
        }}
      
    }}
  
  myYMed<-sqrt((c(DatSumMod$Exposure_Current[length(DatSumMod[,1])],
                  DatSumMod$Exposure_2050_Q90[length(DatSumMod[,1])],
                  DatSumMod$Exposure_2090_Q90[length(DatSumMod[,1])]))+1)
  
  points(myX,myYMed, col="black",bg="black",     pch=21, lwd= 2)
  lines (myX,myYMed, col="black",                lwd= 2, lty=1)
  
  
  # c.ii. percent of current
  
  #plot empty base plot
  myYMed<-c(DatSumMod$Exposure_Current[length(DatSumMod[,1])]/DatSumMod$Exposure_Current[length(DatSumMod[,1])],
            DatSumMod$Exposure_2050_Q90[length(DatSumMod[,1])]/DatSumMod$Exposure_Current[length(DatSumMod[,1])],
            DatSumMod$Exposure_2090_Q90[length(DatSumMod[,1])]/DatSumMod$Exposure_Current[length(DatSumMod[,1])])
  ymax<-max(max(DatSumMod$Exposure_Current/DatSumMod$Exposure_Current,na.rm=TRUE),
            max(DatSumMod$Exposure_2050_Q90/DatSumMod$Exposure_Current,na.rm=TRUE),
            max(DatSumMod$Exposure_2090_Q90/DatSumMod$Exposure_Current,na.rm=TRUE),na.rm=TRUE)
  plot(myX,myYMed, ylim=c(0,ymax),xlim=c(1970,2100),col="white",bg="white", type="b", 
       bty="n", axes=TRUE, pch=21, lwd= 2, lty=1, main="Q90 Exposure [% of Current]")
  
  #plot each species' data
  for (sp in spp){
    n<-which(DatSumMod$Accepted_Species==paste(gsub("_"," ",sp)))
    if(DatSumMod$Exposure_Current[n]>0 & !is.na(DatSumMod$Exposure_Current[n])){
      
      myYMed<-c(DatSumMod$Exposure_Current[n]/DatSumMod$Exposure_Current[n],
                DatSumMod$Exposure_2050_Q90[n]/DatSumMod$Exposure_Current[n],
                DatSumMod$Exposure_2090_Q90[n]/DatSumMod$Exposure_Current[n])
      
      points(myX,myYMed, col = "black", bg="darkgrey",  pch=21, lwd= 1)
      lines (myX,myYMed, col="darkgrey", lwd= 1, lty=1)
      
      if(!is.na(DatSumMod$Exposure_2050_Q90[n])){
        if(DatSumMod$Exposure_2050_Q90[n]>DatSumMod$Exposure_Current[n]){
          lines (myX,myYMed, col="pink", lwd= 1, lty=1)
          print(paste(sp, " increase by 2050"))
        }}
      
      if(!is.na(DatSumMod$Exposure_2090_Q90[n])){
        if(DatSumMod$Exposure_2090_Q90[n]>DatSumMod$Exposure_Current[n]){
          lines (myX,myYMed, col="red", lwd= 1, lty=1)
          print(paste(sp, " increase by 2090"))
        }}
      
    }}
  
  myYMed<-c(DatSumMod$Exposure_Current[length(DatSumMod[,1])]/DatSumMod$Exposure_Current[length(DatSumMod[,1])],
            DatSumMod$Exposure_2050_Q90[length(DatSumMod[,1])]/DatSumMod$Exposure_Current[length(DatSumMod[,1])],
            DatSumMod$Exposure_2090_Q90[length(DatSumMod[,1])]/DatSumMod$Exposure_Current[length(DatSumMod[,1])])
  
  points(myX,myYMed, col="black",bg="black",     pch=21, lwd= 2)
  lines (myX,myYMed, col="black",                lwd= 2, lty=1)
  
  dev.off()
  
}


#####################
# 4. Risk ###########
#####################


{
  png(filename=paste(PlotDir,"/Risk_CCTrends.png",sep=""), units="in", width=10, height=18, res=500) #png smaller file size, A3 for better resolution
  par(mfrow = c(3,2), mar = c(1,1,3,1), oma=c(1,1,5,1), mgp=c(1.2,0.5,0))
  #mar=c(bottom, left, top, right); oma= outer margins; mpg = axis label locations c(3, 1, 0)
  
  myX<-c(1985,2050,2090)
  DatSumMod$Risk_Current[DatSumMod$Risk_Current==0]<-0.01 #avoid dividing by 0 for % change graphs
  
  # a. Median Risk:
  # a.i. raw
  
  #plot empty base plot
  myYMed<- sqrt((   c(DatSumMod$Risk_Current[length(DatSumMod[,1])],
                      DatSumMod$Risk_2050_Median[length(DatSumMod[,1])],
                      DatSumMod$Risk_2090_Median[length(DatSumMod[,1])])   )+1)
  ymax<- sqrt(( max(max(DatSumMod$Risk_Current,na.rm=TRUE),
                    max(DatSumMod$Risk_2050_Median,na.rm=TRUE),
                    max(DatSumMod$Risk_2090_Median,na.rm=TRUE),na.rm=TRUE)     )+1)
  plot(myX,myYMed, ylim=c(0,ymax),xlim=c(1970,2100),col="white",bg="white", type="b", 
       bty="n", axes=TRUE, pch=21, lwd= 2, lty=1, main="Median Risk [sqrt(cell sum)]")
  
  #plot each species' data
  for (sp in spp){
    n<-which(DatSumMod$Accepted_Species==paste(gsub("_"," ",sp)))
    if(DatSumMod$Risk_Current[n]>0 & !is.na(DatSumMod$Risk_Current[n])){
      
      myYMed<-sqrt((  c(DatSumMod$Risk_Current[n],
                        DatSumMod$Risk_2050_Median[n],
                        DatSumMod$Risk_2090_Median[n])     )+1)
      
      points(myX,myYMed, col = "black", bg="darkgrey",  pch=21, lwd= 1)
      lines (myX,myYMed, col="darkgrey", lwd= 1, lty=1)
      
      if(!is.na(DatSumMod$Risk_2050_Median[n])){
        if(DatSumMod$Risk_2050_Median[n]>DatSumMod$Risk_Current[n]){
          lines (myX,myYMed, col="pink", lwd= 1, lty=1)
          print(paste(sp, " increase by 2050"))
        }}
      
      if(!is.na(DatSumMod$Risk_2090_Median[n])){
        if(DatSumMod$Risk_2090_Median[n]>DatSumMod$Risk_Current[n]){
          lines (myX,myYMed, col="red", lwd= 1, lty=1)
          print(paste(sp, " increase by 2090"))
        }}
      
    }}
  
  myYMed<-sqrt((   c(DatSumMod$Risk_Current[length(DatSumMod[,1])],
                     DatSumMod$Risk_2050_Median[length(DatSumMod[,1])],
                     DatSumMod$Risk_2090_Median[length(DatSumMod[,1])])    )+1)
  
  points(myX,myYMed, col="black",bg="black",     pch=21, lwd= 2)
  lines (myX,myYMed, col="black",                lwd= 2, lty=1)
  
  
  # a.ii. percent of current
  
  #plot empty base plot
  myYMed<-c(DatSumMod$Risk_Current[length(DatSumMod[,1])]/DatSumMod$Risk_Current[length(DatSumMod[,1])],
            DatSumMod$Risk_2050_Median[length(DatSumMod[,1])]/DatSumMod$Risk_Current[length(DatSumMod[,1])],
            DatSumMod$Risk_2090_Median[length(DatSumMod[,1])]/DatSumMod$Risk_Current[length(DatSumMod[,1])])
  ymax<-max(max(DatSumMod$Risk_Current/DatSumMod$Risk_Current,na.rm=TRUE),
            max(DatSumMod$Risk_2050_Median/DatSumMod$Risk_Current,na.rm=TRUE),
            max(DatSumMod$Risk_2090_Median/DatSumMod$Risk_Current,na.rm=TRUE),na.rm=TRUE)
  ymax<-ifelse(ymax>10,10,ymax)
  plot(myX,myYMed, ylim=c(0,ymax),xlim=c(1970,2100),col="white",bg="white", type="b", 
       bty="n", axes=TRUE, pch=21, lwd= 2, lty=1, main="Median Risk [% of Current]")
  
  #plot each species' data
  for (sp in spp){
    n<-which(DatSumMod$Accepted_Species==paste(gsub("_"," ",sp)))
    if(DatSumMod$Risk_Current[n]>0 & !is.na(DatSumMod$Risk_Current[n])){
      
      myYMed<-c(DatSumMod$Risk_Current[n]/DatSumMod$Risk_Current[n],
                DatSumMod$Risk_2050_Median[n]/DatSumMod$Risk_Current[n],
                DatSumMod$Risk_2090_Median[n]/DatSumMod$Risk_Current[n])
      
      points(myX,myYMed, col = "black", bg="darkgrey",  pch=21, lwd= 1)
      lines (myX,myYMed, col="darkgrey", lwd= 1, lty=1)
      
      if(!is.na(DatSumMod$Risk_2050_Median[n])){
        if(DatSumMod$Risk_2050_Median[n]>DatSumMod$Risk_Current[n]){
          lines (myX,myYMed, col="pink", lwd= 1, lty=1)
          print(paste(sp, " increase by 2050"))
        }}
      
      if(!is.na(DatSumMod$Risk_2090_Median[n])){
        if(DatSumMod$Risk_2090_Median[n]>DatSumMod$Risk_Current[n]){
          lines (myX,myYMed, col="red", lwd= 1, lty=1)
          print(paste(sp, " increase by 2090"))
        }}
      
    }}
  
  myYMed<-c(DatSumMod$Risk_Current[length(DatSumMod[,1])]/DatSumMod$Risk_Current[length(DatSumMod[,1])],
            DatSumMod$Risk_2050_Median[length(DatSumMod[,1])]/DatSumMod$Risk_Current[length(DatSumMod[,1])],
            DatSumMod$Risk_2090_Median[length(DatSumMod[,1])]/DatSumMod$Risk_Current[length(DatSumMod[,1])])
  
  points(myX,myYMed, col="black",bg="black",     pch=21, lwd= 2)
  lines (myX,myYMed, col="black",                lwd= 2, lty=1)
  
  
  # b. Q10 Risk:
  # b.i. raw
  
  #plot empty base plot
  myYMed<-sqrt((c(DatSumMod$Risk_Current[length(DatSumMod[,1])],
                  DatSumMod$Risk_2050_Q10[length(DatSumMod[,1])],
                  DatSumMod$Risk_2090_Q10[length(DatSumMod[,1])]) )+1)
  ymax<-sqrt((max(max(DatSumMod$Risk_Current,na.rm=TRUE),
                  max(DatSumMod$Risk_2050_Q10,na.rm=TRUE),
                  max(DatSumMod$Risk_2090_Q10,na.rm=TRUE),na.rm=TRUE))+1)
  plot(myX,myYMed, ylim=c(0,ymax),xlim=c(1970,2100),col="white",bg="white", type="b", 
       bty="n", axes=TRUE, pch=21, lwd= 2, lty=1, main="Q10 Risk [sqrt(cell sum)]")
  
  #plot each species' data
  for (sp in spp){
    n<-which(DatSumMod$Accepted_Species==paste(gsub("_"," ",sp)))
    if(DatSumMod$Risk_Current[n]>0 & !is.na(DatSumMod$Risk_Current[n])){
      
      myYMed<-sqrt((c(DatSumMod$Risk_Current[n],
                      DatSumMod$Risk_2050_Q10[n],
                      DatSumMod$Risk_2090_Q10[n]))+1)
      
      points(myX,myYMed, col = "black", bg="darkgrey",  pch=21, lwd= 1)
      lines (myX,myYMed, col="darkgrey", lwd= 1, lty=1)
      
      if(!is.na(DatSumMod$Risk_2050_Q10[n])){
        if(DatSumMod$Risk_2050_Q10[n]>DatSumMod$Risk_Current[n]){
          lines (myX,myYMed, col="pink", lwd= 1, lty=1)
          print(paste(sp, " increase by 2050"))
        }}
      
      if(!is.na(DatSumMod$Risk_2090_Q10[n])){
        if(DatSumMod$Risk_2090_Q10[n]>DatSumMod$Risk_Current[n]){
          lines (myX,myYMed, col="red", lwd= 1, lty=1)
          print(paste(sp, " increase by 2090"))
        }}
      
    }}
  
  myYMed<-sqrt((c(DatSumMod$Risk_Current[length(DatSumMod[,1])],
                  DatSumMod$Risk_2050_Q10[length(DatSumMod[,1])],
                  DatSumMod$Risk_2090_Q10[length(DatSumMod[,1])]))+1)
  
  points(myX,myYMed, col="black",bg="black",     pch=21, lwd= 2)
  lines (myX,myYMed, col="black",                lwd= 2, lty=1)
  
  
  # b.ii. percent of current
  
  #plot empty base plot
  myYMed<-c(DatSumMod$Risk_Current[length(DatSumMod[,1])]/DatSumMod$Risk_Current[length(DatSumMod[,1])],
            DatSumMod$Risk_2050_Q10[length(DatSumMod[,1])]/DatSumMod$Risk_Current[length(DatSumMod[,1])],
            DatSumMod$Risk_2090_Q10[length(DatSumMod[,1])]/DatSumMod$Risk_Current[length(DatSumMod[,1])])
  ymax<-max(max(DatSumMod$Risk_Current/DatSumMod$Risk_Current,na.rm=TRUE),
            max(DatSumMod$Risk_2050_Q10/DatSumMod$Risk_Current,na.rm=TRUE),
            max(DatSumMod$Risk_2090_Q10/DatSumMod$Risk_Current,na.rm=TRUE),na.rm=TRUE)
  ymax<-ifelse(ymax>10,10,ymax)
  plot(myX,myYMed, ylim=c(0,ymax),xlim=c(1970,2100),col="white",bg="white", type="b", 
       bty="n", axes=TRUE, pch=21, lwd= 2, lty=1, main="Q10 Risk [% of Current]")
  
  #plot each species' data
  for (sp in spp){
    n<-which(DatSumMod$Accepted_Species==paste(gsub("_"," ",sp)))
    if(DatSumMod$Risk_Current[n]>0 & !is.na(DatSumMod$Risk_Current[n])){
      
      myYMed<-c(DatSumMod$Risk_Current[n]/DatSumMod$Risk_Current[n],
                DatSumMod$Risk_2050_Q10[n]/DatSumMod$Risk_Current[n],
                DatSumMod$Risk_2090_Q10[n]/DatSumMod$Risk_Current[n])
      
      points(myX,myYMed, col = "black", bg="darkgrey",  pch=21, lwd= 1)
      lines (myX,myYMed, col="darkgrey", lwd= 1, lty=1)
      
      if(!is.na(DatSumMod$Risk_2050_Q10[n])){
        if(DatSumMod$Risk_2050_Q10[n]>DatSumMod$Risk_Current[n]){
          lines (myX,myYMed, col="pink", lwd= 1, lty=1)
          print(paste(sp, " increase by 2050"))
        }}
      
      if(!is.na(DatSumMod$Risk_2090_Q10[n])){
        if(DatSumMod$Risk_2090_Q10[n]>DatSumMod$Risk_Current[n]){
          lines (myX,myYMed, col="red", lwd= 1, lty=1)
          print(paste(sp, " increase by 2090"))
          
        }}
      
    }}
  
  myYMed<-c(DatSumMod$Risk_Current[length(DatSumMod[,1])]/DatSumMod$Risk_Current[length(DatSumMod[,1])],
            DatSumMod$Risk_2050_Q10[length(DatSumMod[,1])]/DatSumMod$Risk_Current[length(DatSumMod[,1])],
            DatSumMod$Risk_2090_Q10[length(DatSumMod[,1])]/DatSumMod$Risk_Current[length(DatSumMod[,1])])
  
  points(myX,myYMed, col="black",bg="black",     pch=21, lwd= 2)
  lines (myX,myYMed, col="black",                lwd= 2, lty=1)
  
  # c. Q90 Risk:
  # c.i. raw
  
  #plot empty base plot
  myYMed<-sqrt((c(DatSumMod$Risk_Current[length(DatSumMod[,1])],
                  DatSumMod$Risk_2050_Q90[length(DatSumMod[,1])],
                  DatSumMod$Risk_2090_Q90[length(DatSumMod[,1])]))+1)
  ymax<-sqrt((max(max(DatSumMod$Risk_Current,na.rm=TRUE),
                  max(DatSumMod$Risk_2050_Q90,na.rm=TRUE),
                  max(DatSumMod$Risk_2090_Q90,na.rm=TRUE),na.rm=TRUE))+1)
  plot(myX,myYMed, ylim=c(0,ymax),xlim=c(1970,2100),col="white",bg="white", type="b", 
       bty="n", axes=TRUE, pch=21, lwd= 2, lty=1, main="Q90 Risk [sqrt(cell sum)]")
  
  #plot each species' data
  for (sp in spp){
    n<-which(DatSumMod$Accepted_Species==paste(gsub("_"," ",sp)))
    if(DatSumMod$Risk_Current[n]>0 & !is.na(DatSumMod$Risk_Current[n])){
      
      myYMed<-sqrt((c(DatSumMod$Risk_Current[n],
                      DatSumMod$Risk_2050_Q90[n],
                      DatSumMod$Risk_2090_Q90[n]))+1)
      
      points(myX,myYMed, col = "black", bg="darkgrey",  pch=21, lwd= 1)
      lines (myX,myYMed, col="darkgrey", lwd= 1, lty=1)
      
      if(!is.na(DatSumMod$Risk_2050_Q90[n])){
        if(DatSumMod$Risk_2050_Q90[n]>DatSumMod$Risk_Current[n]){
          lines (myX,myYMed, col="pink", lwd= 1, lty=1)
          print(paste(sp, " increase by 2050"))
        }}
      
      if(!is.na(DatSumMod$Risk_2090_Q90[n])){
        if(DatSumMod$Risk_2090_Q90[n]>DatSumMod$Risk_Current[n]){
          lines (myX,myYMed, col="red", lwd= 1, lty=1)
          print(paste(sp, " increase by 2090"))
        }}
      
    }}
  
  myYMed<-sqrt((c(DatSumMod$Risk_Current[length(DatSumMod[,1])],
                  DatSumMod$Risk_2050_Q90[length(DatSumMod[,1])],
                  DatSumMod$Risk_2090_Q90[length(DatSumMod[,1])]))+1)
  
  points(myX,myYMed, col="black",bg="black",     pch=21, lwd= 2)
  lines (myX,myYMed, col="black",                lwd= 2, lty=1)
  
  
  # c.ii. percent of current
  
  #plot empty base plot
  myYMed<-c(DatSumMod$Risk_Current[length(DatSumMod[,1])]/DatSumMod$Risk_Current[length(DatSumMod[,1])],
            DatSumMod$Risk_2050_Q90[length(DatSumMod[,1])]/DatSumMod$Risk_Current[length(DatSumMod[,1])],
            DatSumMod$Risk_2090_Q90[length(DatSumMod[,1])]/DatSumMod$Risk_Current[length(DatSumMod[,1])])
  ymax<-max(max(DatSumMod$Risk_Current/DatSumMod$Risk_Current,na.rm=TRUE),
            max(DatSumMod$Risk_2050_Q90/DatSumMod$Risk_Current,na.rm=TRUE),
            max(DatSumMod$Risk_2090_Q90/DatSumMod$Risk_Current,na.rm=TRUE),na.rm=TRUE)
 ymax<-ifelse(ymax>10,10,ymax)
   plot(myX,myYMed, ylim=c(0,ymax),xlim=c(1970,2100),col="white",bg="white", type="b", 
       bty="n", axes=TRUE, pch=21, lwd= 2, lty=1, main="Q90 Risk [% of Current]")
  
  #plot each species' data
  for (sp in spp){
    n<-which(DatSumMod$Accepted_Species==paste(gsub("_"," ",sp)))
    if(DatSumMod$Risk_Current[n]>0 & !is.na(DatSumMod$Risk_Current[n])){
      
      myYMed<-c(DatSumMod$Risk_Current[n]/DatSumMod$Risk_Current[n],
                DatSumMod$Risk_2050_Q90[n]/DatSumMod$Risk_Current[n],
                DatSumMod$Risk_2090_Q90[n]/DatSumMod$Risk_Current[n])
      
      points(myX,myYMed, col = "black", bg="darkgrey",  pch=21, lwd= 1)
      lines (myX,myYMed, col="darkgrey", lwd= 1, lty=1)
      
      if(!is.na(DatSumMod$Risk_2050_Q90[n])){
        if(DatSumMod$Risk_2050_Q90[n]>DatSumMod$Risk_Current[n]){
          lines (myX,myYMed, col="pink", lwd= 1, lty=1)
          print(paste(sp, " increase by 2050"))
        }}
      
      if(!is.na(DatSumMod$Risk_2090_Q90[n])){
        if(DatSumMod$Risk_2090_Q90[n]>DatSumMod$Risk_Current[n]){
          lines (myX,myYMed, col="red", lwd= 1, lty=1)
          print(paste(sp, " increase by 2090"))
        }}
      
    }}
  
  myYMed<-c(DatSumMod$Risk_Current[length(DatSumMod[,1])]/DatSumMod$Risk_Current[length(DatSumMod[,1])],
            DatSumMod$Risk_2050_Q90[length(DatSumMod[,1])]/DatSumMod$Risk_Current[length(DatSumMod[,1])],
            DatSumMod$Risk_2090_Q90[length(DatSumMod[,1])]/DatSumMod$Risk_Current[length(DatSumMod[,1])])
  
  points(myX,myYMed, col="black",bg="black",     pch=21, lwd= 2)
  lines (myX,myYMed, col="black",                lwd= 2, lty=1)
  
  dev.off()
  
}


#####################
# 5. Impact #########
#####################


{
  png(filename=paste(PlotDir,"/Impact_CCTrends.png",sep=""), units="in", width=10, height=18, res=500) #png smaller file size, A3 for better resolution
  par(mfrow = c(3,2), mar = c(1,1,3,1), oma=c(1,1,5,1), mgp=c(1.2,0.5,0))
  #mar=c(bottom, left, top, right); oma= outer margins; mpg = axis label locations c(3, 1, 0)
  
  myX<-c(1985,2050,2090)
  DatSumMod$Impact_Current[DatSumMod$Impact_Current==0]<-0.01 #avoid dividing by 0 for % change graphs
  
  # a. Median Impact:
  # a.i. raw
  
  #plot empty base plot
  myYMed<- sqrt((   c(DatSumMod$Impact_Current[length(DatSumMod[,1])],
                      DatSumMod$Impact_2050_Median[length(DatSumMod[,1])],
                      DatSumMod$Impact_2090_Median[length(DatSumMod[,1])])   )+1)
  ymax<- sqrt(( max(max(DatSumMod$Impact_Current,na.rm=TRUE),
                    max(DatSumMod$Impact_2050_Median,na.rm=TRUE),
                    max(DatSumMod$Impact_2090_Median,na.rm=TRUE),na.rm=TRUE)     )+1)
  plot(myX,myYMed, ylim=c(0,ymax),xlim=c(1970,2100),col="white",bg="white", type="b", 
       bty="n", axes=TRUE, pch=21, lwd= 2, lty=1, main="Median Impact [sqrt(cell sum)]")
  
  #plot each species' data
  for (sp in spp){
    n<-which(DatSumMod$Accepted_Species==paste(gsub("_"," ",sp)))
    if(DatSumMod$Impact_Current[n]>0 & !is.na(DatSumMod$Impact_Current[n])){
      
      myYMed<-sqrt((  c(DatSumMod$Impact_Current[n],
                        DatSumMod$Impact_2050_Median[n],
                        DatSumMod$Impact_2090_Median[n])     )+1)
      
      points(myX,myYMed, col = "black", bg="darkgrey",  pch=21, lwd= 1)
      lines (myX,myYMed, col="darkgrey", lwd= 1, lty=1)
      
      if(!is.na(DatSumMod$Impact_2050_Median[n])){
        if(DatSumMod$Impact_2050_Median[n]>DatSumMod$Impact_Current[n]){
          lines (myX,myYMed, col="pink", lwd= 1, lty=1)
          print(paste(sp, " increase by 2050"))
        }}
      
      if(!is.na(DatSumMod$Impact_2090_Median[n])){
        if(DatSumMod$Impact_2090_Median[n]>DatSumMod$Impact_Current[n]){
          lines (myX,myYMed, col="red", lwd= 1, lty=1)
          print(paste(sp, " increase by 2090"))
        }}
      
    }}
  
  myYMed<-sqrt((   c(DatSumMod$Impact_Current[length(DatSumMod[,1])],
                     DatSumMod$Impact_2050_Median[length(DatSumMod[,1])],
                     DatSumMod$Impact_2090_Median[length(DatSumMod[,1])])    )+1)
  
  points(myX,myYMed, col="black",bg="black",     pch=21, lwd= 2)
  lines (myX,myYMed, col="black",                lwd= 2, lty=1)
  
  
  # a.ii. percent of current
  
  #plot empty base plot
  myYMed<-c(DatSumMod$Impact_Current[length(DatSumMod[,1])]/DatSumMod$Impact_Current[length(DatSumMod[,1])],
            DatSumMod$Impact_2050_Median[length(DatSumMod[,1])]/DatSumMod$Impact_Current[length(DatSumMod[,1])],
            DatSumMod$Impact_2090_Median[length(DatSumMod[,1])]/DatSumMod$Impact_Current[length(DatSumMod[,1])])
  ymax<-max(max(DatSumMod$Impact_Current/DatSumMod$Impact_Current,na.rm=TRUE),
            max(DatSumMod$Impact_2050_Median/DatSumMod$Impact_Current,na.rm=TRUE),
            max(DatSumMod$Impact_2090_Median/DatSumMod$Impact_Current,na.rm=TRUE),na.rm=TRUE)
  ymax<-ifelse(ymax>10,10,ymax)
  plot(myX,myYMed, ylim=c(0,ymax),xlim=c(1970,2100),col="white",bg="white", type="b", 
       bty="n", axes=TRUE, pch=21, lwd= 2, lty=1, main="Median Impact [% of Current]")
  
  #plot each species' data
  for (sp in spp){
    n<-which(DatSumMod$Accepted_Species==paste(gsub("_"," ",sp)))
    if(DatSumMod$Impact_Current[n]>0 & !is.na(DatSumMod$Impact_Current[n])){
      
      myYMed<-c(DatSumMod$Impact_Current[n]/DatSumMod$Impact_Current[n],
                DatSumMod$Impact_2050_Median[n]/DatSumMod$Impact_Current[n],
                DatSumMod$Impact_2090_Median[n]/DatSumMod$Impact_Current[n])
      
      points(myX,myYMed, col = "black", bg="darkgrey",  pch=21, lwd= 1)
      lines (myX,myYMed, col="darkgrey", lwd= 1, lty=1)
      
      if(!is.na(DatSumMod$Impact_2050_Median[n])){
        if(DatSumMod$Impact_2050_Median[n]>DatSumMod$Impact_Current[n]){
          lines (myX,myYMed, col="pink", lwd= 1, lty=1)
          print(paste(sp, " increase by 2050"))
        }}
      
      if(!is.na(DatSumMod$Impact_2090_Median[n])){
        if(DatSumMod$Impact_2090_Median[n]>DatSumMod$Impact_Current[n]){
          lines (myX,myYMed, col="red", lwd= 1, lty=1)
          print(paste(sp, " increase by 2090"))
        }}
      
    }}
  
  myYMed<-c(DatSumMod$Impact_Current[length(DatSumMod[,1])]/DatSumMod$Impact_Current[length(DatSumMod[,1])],
            DatSumMod$Impact_2050_Median[length(DatSumMod[,1])]/DatSumMod$Impact_Current[length(DatSumMod[,1])],
            DatSumMod$Impact_2090_Median[length(DatSumMod[,1])]/DatSumMod$Impact_Current[length(DatSumMod[,1])])
  
  points(myX,myYMed, col="black",bg="black",     pch=21, lwd= 2)
  lines (myX,myYMed, col="black",                lwd= 2, lty=1)
  
  
  # b. Q10 Impact:
  # b.i. raw
  
  #plot empty base plot
  myYMed<-sqrt((c(DatSumMod$Impact_Current[length(DatSumMod[,1])],
                  DatSumMod$Impact_2050_Q10[length(DatSumMod[,1])],
                  DatSumMod$Impact_2090_Q10[length(DatSumMod[,1])]) )+1)
  ymax<-sqrt((max(max(DatSumMod$Impact_Current,na.rm=TRUE),
                  max(DatSumMod$Impact_2050_Q10,na.rm=TRUE),
                  max(DatSumMod$Impact_2090_Q10,na.rm=TRUE),na.rm=TRUE))+1)
  plot(myX,myYMed, ylim=c(0,ymax),xlim=c(1970,2100),col="white",bg="white", type="b", 
       bty="n", axes=TRUE, pch=21, lwd= 2, lty=1, main="Q10 Impact [sqrt(cell sum)]")
  
  #plot each species' data
  for (sp in spp){
    n<-which(DatSumMod$Accepted_Species==paste(gsub("_"," ",sp)))
    if(DatSumMod$Impact_Current[n]>0 & !is.na(DatSumMod$Impact_Current[n])){
      
      myYMed<-sqrt((c(DatSumMod$Impact_Current[n],
                      DatSumMod$Impact_2050_Q10[n],
                      DatSumMod$Impact_2090_Q10[n]))+1)
      
      points(myX,myYMed, col = "black", bg="darkgrey",  pch=21, lwd= 1)
      lines (myX,myYMed, col="darkgrey", lwd= 1, lty=1)
      
      if(!is.na(DatSumMod$Impact_2050_Q10[n])){
        if(DatSumMod$Impact_2050_Q10[n]>DatSumMod$Impact_Current[n]){
          lines (myX,myYMed, col="pink", lwd= 1, lty=1)
          print(paste(sp, " increase by 2050"))
        }}
      
      if(!is.na(DatSumMod$Impact_2090_Q10[n])){
        if(DatSumMod$Impact_2090_Q10[n]>DatSumMod$Impact_Current[n]){
          lines (myX,myYMed, col="red", lwd= 1, lty=1)
          print(paste(sp, " increase by 2090"))
        }}
      
    }}
  
  myYMed<-sqrt((c(DatSumMod$Impact_Current[length(DatSumMod[,1])],
                  DatSumMod$Impact_2050_Q10[length(DatSumMod[,1])],
                  DatSumMod$Impact_2090_Q10[length(DatSumMod[,1])]))+1)
  
  points(myX,myYMed, col="black",bg="black",     pch=21, lwd= 2)
  lines (myX,myYMed, col="black",                lwd= 2, lty=1)
  
  
  # b.ii. percent of current
  
  #plot empty base plot
  myYMed<-c(DatSumMod$Impact_Current[length(DatSumMod[,1])]/DatSumMod$Impact_Current[length(DatSumMod[,1])],
            DatSumMod$Impact_2050_Q10[length(DatSumMod[,1])]/DatSumMod$Impact_Current[length(DatSumMod[,1])],
            DatSumMod$Impact_2090_Q10[length(DatSumMod[,1])]/DatSumMod$Impact_Current[length(DatSumMod[,1])])
  ymax<-max(max(DatSumMod$Impact_Current/DatSumMod$Impact_Current,na.rm=TRUE),
            max(DatSumMod$Impact_2050_Q10/DatSumMod$Impact_Current,na.rm=TRUE),
            max(DatSumMod$Impact_2090_Q10/DatSumMod$Impact_Current,na.rm=TRUE),na.rm=TRUE)
  ymax<-ifelse(ymax>10,10,ymax)
  plot(myX,myYMed, ylim=c(0,ymax),xlim=c(1970,2100),col="white",bg="white", type="b", 
       bty="n", axes=TRUE, pch=21, lwd= 2, lty=1, main="Q10 Impact [% of Current]")
  
  #plot each species' data
  for (sp in spp){
    n<-which(DatSumMod$Accepted_Species==paste(gsub("_"," ",sp)))
    if(DatSumMod$Impact_Current[n]>0 & !is.na(DatSumMod$Impact_Current[n])){
      
      myYMed<-c(DatSumMod$Impact_Current[n]/DatSumMod$Impact_Current[n],
                DatSumMod$Impact_2050_Q10[n]/DatSumMod$Impact_Current[n],
                DatSumMod$Impact_2090_Q10[n]/DatSumMod$Impact_Current[n])
      
      points(myX,myYMed, col = "black", bg="darkgrey",  pch=21, lwd= 1)
      lines (myX,myYMed, col="darkgrey", lwd= 1, lty=1)
      
      if(!is.na(DatSumMod$Impact_2050_Q10[n])){
        if(DatSumMod$Impact_2050_Q10[n]>DatSumMod$Impact_Current[n]){
          lines (myX,myYMed, col="pink", lwd= 1, lty=1)
          print(paste(sp, " increase by 2050"))
        }}
      
      if(!is.na(DatSumMod$Impact_2090_Q10[n])){
        if(DatSumMod$Impact_2090_Q10[n]>DatSumMod$Impact_Current[n]){
          lines (myX,myYMed, col="red", lwd= 1, lty=1)
          print(paste(sp, " increase by 2090"))
          
        }}
      
    }}
  
  myYMed<-c(DatSumMod$Impact_Current[length(DatSumMod[,1])]/DatSumMod$Impact_Current[length(DatSumMod[,1])],
            DatSumMod$Impact_2050_Q10[length(DatSumMod[,1])]/DatSumMod$Impact_Current[length(DatSumMod[,1])],
            DatSumMod$Impact_2090_Q10[length(DatSumMod[,1])]/DatSumMod$Impact_Current[length(DatSumMod[,1])])
  
  points(myX,myYMed, col="black",bg="black",     pch=21, lwd= 2)
  lines (myX,myYMed, col="black",                lwd= 2, lty=1)
  
  # c. Q90 Impact:
  # c.i. raw
  
  #plot empty base plot
  myYMed<-sqrt((c(DatSumMod$Impact_Current[length(DatSumMod[,1])],
                  DatSumMod$Impact_2050_Q90[length(DatSumMod[,1])],
                  DatSumMod$Impact_2090_Q90[length(DatSumMod[,1])]))+1)
  ymax<-sqrt((max(max(DatSumMod$Impact_Current,na.rm=TRUE),
                  max(DatSumMod$Impact_2050_Q90,na.rm=TRUE),
                  max(DatSumMod$Impact_2090_Q90,na.rm=TRUE),na.rm=TRUE))+1)
  plot(myX,myYMed, ylim=c(0,ymax),xlim=c(1970,2100),col="white",bg="white", type="b", 
       bty="n", axes=TRUE, pch=21, lwd= 2, lty=1, main="Q90 Impact [sqrt(cell sum)]")
  
  #plot each species' data
  for (sp in spp){
    n<-which(DatSumMod$Accepted_Species==paste(gsub("_"," ",sp)))
    if(DatSumMod$Impact_Current[n]>0 & !is.na(DatSumMod$Impact_Current[n])){
      
      myYMed<-sqrt((c(DatSumMod$Impact_Current[n],
                      DatSumMod$Impact_2050_Q90[n],
                      DatSumMod$Impact_2090_Q90[n]))+1)
      
      points(myX,myYMed, col = "black", bg="darkgrey",  pch=21, lwd= 1)
      lines (myX,myYMed, col="darkgrey", lwd= 1, lty=1)
      
      if(!is.na(DatSumMod$Impact_2050_Q90[n])){
        if(DatSumMod$Impact_2050_Q90[n]>DatSumMod$Impact_Current[n]){
          lines (myX,myYMed, col="pink", lwd= 1, lty=1)
          print(paste(sp, " increase by 2050"))
        }}
      
      if(!is.na(DatSumMod$Impact_2090_Q90[n])){
        if(DatSumMod$Impact_2090_Q90[n]>DatSumMod$Impact_Current[n]){
          lines (myX,myYMed, col="red", lwd= 1, lty=1)
          print(paste(sp, " increase by 2090"))
        }}
      
    }}
  
  myYMed<-sqrt((c(DatSumMod$Impact_Current[length(DatSumMod[,1])],
                  DatSumMod$Impact_2050_Q90[length(DatSumMod[,1])],
                  DatSumMod$Impact_2090_Q90[length(DatSumMod[,1])]))+1)
  
  points(myX,myYMed, col="black",bg="black",     pch=21, lwd= 2)
  lines (myX,myYMed, col="black",                lwd= 2, lty=1)
  
  
  # c.ii. percent of current
  
  #plot empty base plot
  myYMed<-c(DatSumMod$Impact_Current[length(DatSumMod[,1])]/DatSumMod$Impact_Current[length(DatSumMod[,1])],
            DatSumMod$Impact_2050_Q90[length(DatSumMod[,1])]/DatSumMod$Impact_Current[length(DatSumMod[,1])],
            DatSumMod$Impact_2090_Q90[length(DatSumMod[,1])]/DatSumMod$Impact_Current[length(DatSumMod[,1])])
  ymax<-max(max(DatSumMod$Impact_Current/DatSumMod$Impact_Current,na.rm=TRUE),
            max(DatSumMod$Impact_2050_Q90/DatSumMod$Impact_Current,na.rm=TRUE),
            max(DatSumMod$Impact_2090_Q90/DatSumMod$Impact_Current,na.rm=TRUE),na.rm=TRUE)
  ymax<-ifelse(ymax>10,10,ymax)
  plot(myX,myYMed, ylim=c(0,ymax),xlim=c(1970,2100),col="white",bg="white", type="b", 
       bty="n", axes=TRUE, pch=21, lwd= 2, lty=1, main="Q90 Impact [% of Current]")
  
  #plot each species' data
  for (sp in spp){
    n<-which(DatSumMod$Accepted_Species==paste(gsub("_"," ",sp)))
    if(DatSumMod$Impact_Current[n]>0 & !is.na(DatSumMod$Impact_Current[n])){
      
      myYMed<-c(DatSumMod$Impact_Current[n]/DatSumMod$Impact_Current[n],
                DatSumMod$Impact_2050_Q90[n]/DatSumMod$Impact_Current[n],
                DatSumMod$Impact_2090_Q90[n]/DatSumMod$Impact_Current[n])
      
      points(myX,myYMed, col = "black", bg="darkgrey",  pch=21, lwd= 1)
      lines (myX,myYMed, col="darkgrey", lwd= 1, lty=1)
      
      if(!is.na(DatSumMod$Impact_2050_Q90[n])){
        if(DatSumMod$Impact_2050_Q90[n]>DatSumMod$Impact_Current[n]){
          lines (myX,myYMed, col="pink", lwd= 1, lty=1)
          print(paste(sp, " increase by 2050"))
        }}
      
      if(!is.na(DatSumMod$Impact_2090_Q90[n])){
        if(DatSumMod$Impact_2090_Q90[n]>DatSumMod$Impact_Current[n]){
          lines (myX,myYMed, col="red", lwd= 1, lty=1)
          print(paste(sp, " increase by 2090"))
        }}
      
    }}
  
  myYMed<-c(DatSumMod$Impact_Current[length(DatSumMod[,1])]/DatSumMod$Impact_Current[length(DatSumMod[,1])],
            DatSumMod$Impact_2050_Q90[length(DatSumMod[,1])]/DatSumMod$Impact_Current[length(DatSumMod[,1])],
            DatSumMod$Impact_2090_Q90[length(DatSumMod[,1])]/DatSumMod$Impact_Current[length(DatSumMod[,1])])
  
  points(myX,myYMed, col="black",bg="black",     pch=21, lwd= 2)
  lines (myX,myYMed, col="black",                lwd= 2, lty=1)
  
  dev.off()
  
}


#####################
# 6. Distance Shift #
#####################



#####################
# 7. Latitude Shift #
#####################

{
png(filename=paste(PlotDir,"/Latitudes.png",sep=""), units="in", width=10, height=10, res=500) #png smaller file size, A3 for better resolution
par(mfrow = c(1,1), mar = c(1,1,3,1), oma=c(1,1,5,1), mgp=c(1.2,0.5,0))

plot(DatSumMod$lat_Current, DatSumMod$lat_2050_Median-DatSumMod$lat_Current, ylim=c(-10,10),xlim=c(-50,50),col="black",bg="white", 
     type="p", bty="n", axes=TRUE, pch=21, lwd= 2, main="Distribution center latitudes",cex=1.5)

points(DatSumMod$lat_Current,DatSumMod$lat_2090_Median-DatSumMod$lat_Current, col="black",     pch=8, lwd= 2)
lines(c(-50,50),c(0,0),lwd=2)
lines(c(0,0),c(-10,10),lwd=2)
dev.off()
}

###############################
# 8a. Circular Direction Shift #
###############################


#########################################
# 8b. Map of direction & distance shift #
#########################################

#plot World Map

#for each species

#create a vector line from current to new range center 2050 & plot
#create a vector line from 2050 to 2090 range center & plot




###############################
# 9. AUC ######################
###############################

{
png(filename=paste(PlotDir,"/AUCs.png",sep=""), units="in", width=10, height=10, res=500) #png smaller file size, A3 for better resolution
par(mfrow = c(1,1), mar = c(1,1,3,1), oma=c(1,1,5,1), mgp=c(1.2,0.5,0))

plot(log(DatSumMod$MUrecs+1), DatSumMod$TrainingAUC, ylim=c(0,1),xlim=c(0,max(log(DatSumMod$MUrecs+1))),col="black",bg="white", 
     type="p", bty="n", axes=TRUE, pch=21, lwd= 2, main="Model AUCs",cex=1.5)

points(log(DatSumMod$MUrecs+1),DatSumMod$TestAUC, col="black",     pch=8, lwd= 2)

dev.off()
}


################################


myYMed<-c(DatSumMod$RangesizeSQKM_Current[length(DatSumMod[,1])]/DatSumMod$RangesizeSQKM_Current[length(DatSumMod[,1])],DatSumMod$RangesizeSQKM_2050_Median[length(DatSumMod[,1])]/DatSumMod$RangesizeSQKM_Current[length(DatSumMod[,1])],DatSumMod$RangesizeSQKM_2090_Median[length(DatSumMod[,1])]/DatSumMod$RangesizeSQKM_Current[length(DatSumMod[,1])])

plot(myX,myYMed, ylim=c(-0.5,3),xlim=c(1970,2100),col="white",bg="white", type="b", bty="n", axes=TRUE, pch=21, lwd= 2, lty=1, main="Range Size [sqkm]")

for (sp in spp){
n<-which(DatSumMod$Accepted_Species==paste(gsub("_"," ",sp)))
if(DatSumMod$RangesizeSQKM_Current[n]>0 & !is.na(DatSumMod$RangesizeSQKM_Current[n])){
  
myYMed<-c(DatSumMod$RangesizeSQKM_Current[n]/DatSumMod$RangesizeSQKM_Current[n],DatSumMod$RangesizeSQKM_2050_Median[n]/DatSumMod$RangesizeSQKM_Current[n],DatSumMod$RangesizeSQKM_2090_Median[n]/DatSumMod$RangesizeSQKM_Current[n])
# myYQ10<-c(DatSumMod$RangesizeSQKM_Current[n]/DatSumMod$RangesizeSQKM_Current[n],DatSumMod$RangesizeSQKM_2050_Q10[n]/DatSumMod$RangesizeSQKM_Current[n],DatSumMod$RangesizeSQKM_2090_Q10[n]/DatSumMod$RangesizeSQKM_Current[n])
# myYQ90<-c(DatSumMod$RangesizeSQKM_Current[n]/DatSumMod$RangesizeSQKM_Current[n],DatSumMod$RangesizeSQKM_2050_Q90[n]/DatSumMod$RangesizeSQKM_Current[n],DatSumMod$RangesizeSQKM_2090_Q90[n]/DatSumMod$RangesizeSQKM_Current[n])

points(myX,myYMed, col = "black", bg="darkgrey",  pch=21, lwd= 1)
lines (myX,myYMed, col="darkgrey", lwd= 1, lty=1)
# points(myX,myYQ10, col = "navy",bg="lightblue1",  pch=21, lwd= 1)
# lines (myX,myYQ10, col="lightblue", lwd= 1, lty=2)
# points(myX,myYQ90, col = "darkred",bg="pink",     pch=21, lwd= 1)
# lines (myX,myYQ90, col="pink", lwd= 1, lty=3)

if(DatSumMod$RangesizeSQKM_2050_Median[n]>DatSumMod$RangesizeSQKM_Current[n]){
  lines (myX,myYMed, col="pink", lwd= 1, lty=1)
  print(paste(sp, " increase by 2050"))
}

if(DatSumMod$RangesizeSQKM_2090_Median[n]>DatSumMod$RangesizeSQKM_Current[n]){
  lines (myX,myYMed, col="red", lwd= 1, lty=1)
  print(paste(sp, " increase by 2090"))
  
}


}}

myYMed<-c(DatSumMod$RangesizeSQKM_Current[length(DatSumMod[,1])]/DatSumMod$RangesizeSQKM_Current[length(DatSumMod[,1])],DatSumMod$RangesizeSQKM_2050_Median[length(DatSumMod[,1])]/DatSumMod$RangesizeSQKM_Current[length(DatSumMod[,1])],DatSumMod$RangesizeSQKM_2090_Median[length(DatSumMod[,1])]/DatSumMod$RangesizeSQKM_Current[length(DatSumMod[,1])])
# myYQ10<-c(DatSumMod$RangesizeSQKM_Current[length(DatSumMod[,1])]/DatSumMod$RangesizeSQKM_Current[length(DatSumMod[,1])],DatSumMod$RangesizeSQKM_2050_Q10[length(DatSumMod[,1])]/DatSumMod$RangesizeSQKM_Current[length(DatSumMod[,1])],   DatSumMod$RangesizeSQKM_2090_Q10[length(DatSumMod[,1])]/DatSumMod$RangesizeSQKM_Current[length(DatSumMod[,1])])
# myYQ90<-c(DatSumMod$RangesizeSQKM_Current[length(DatSumMod[,1])]/DatSumMod$RangesizeSQKM_Current[length(DatSumMod[,1])],DatSumMod$RangesizeSQKM_2050_Q90[length(DatSumMod[,1])]/DatSumMod$RangesizeSQKM_Current[length(DatSumMod[,1])],   DatSumMod$RangesizeSQKM_2090_Q90[length(DatSumMod[,1])]/DatSumMod$RangesizeSQKM_Current[length(DatSumMod[,1])])

points(myX,myYMed, col="black",bg="black",     pch=21, lwd= 2)
lines (myX,myYMed, col="black",                lwd= 2, lty=1)
# points(myX,myYQ10, col="navy",bg="navy",       pch=21, lwd= 2)
# lines (myX,myYQ10, col="navy",                 lwd= 2, lty=2)
# points(myX,myYQ90, col="darkred",bg="darkred", pch=21, lwd= 2)
# lines (myX,myYQ90, col="darkred",              lwd= 2, lty=3)

dev.off()


#########################################
png(filename=paste(PlotDir,"/ClimatChangePlots2.png",sep=""), units="in", width=14, height=14, res=500) #png smaller file size, A3 for better resolution
plot(DatSumMod$RangesizeSQKM_Current,DatSumMod$RangesizeSQKM_2050_Median,col="black",bg="black", xlim=c(0,25000000),ylim=c(0,25000000), type="p", bty="n", axes=TRUE, pch=21, cex=1.2, main="Range Size [sqkm]")
points(DatSumMod$RangesizeSQKM_Current,DatSumMod$RangesizeSQKM_2090_Median,col="grey",bg="grey", pch=21, cex=0.8)
lines(c(0,25000000),c(0,25000000))
dev.off()


#############################################################
#############################################################
#############################################################
#############################################################
#long version of vector calculations

#################################################
# !!!dissolve polygons in ArcGIS first!!!!! #####
#################################################

# #read dissolved range polygons
# SDMshapes_Current <- st_read(paste(IndSppOutDir,"/Current_SDMfeatureAll_dissolved.shp",sep="")) #dissolve in ArcGIS first, easier
# SDMshapes_2050_Median <- st_read(paste(IndSppOutDir,"/2050_Median_SDMfeatureAll_dissolved.shp",sep="")) #dissolve in ArcGIS first, easier
# SDMshapes_2090_Median <- st_read(paste(IndSppOutDir,"/2090_Median_SDMfeatureAll_dissolved.shp",sep="")) #dissolve in ArcGIS first, easier
# # SDMshapes_2050_Q10 <- st_read(paste(IndSppOutDir,"/2050_Q10_SDMfeatureAll_dissolved.shp",sep="")) #dissolve in ArcGIS first, easier
# # SDMshapes_2090_Q10 <- st_read(paste(IndSppOutDir,"/2090_Q10_SDMfeatureAll_dissolved.shp",sep="")) #dissolve in ArcGIS first, easier
# # SDMshapes_2050_Q90 <- st_read(paste(IndSppOutDir,"/2050_Q90_SDMfeatureAll_dissolved.shp",sep="")) #dissolve in ArcGIS first, easier
# # SDMshapes_2090_Q90 <- st_read(paste(IndSppOutDir,"/2090_Q90_SDMfeatureAll_dissolved.shp",sep="")) #dissolve in ArcGIS first, easier
# 
# #calculate SDM distribution centroids:
# cent_Current <- as.data.frame(st_centroid(SDMshapes_Current))
# cent_Current<-cbind(cent_Current,st_coordinates(st_centroid(SDMshapes_Current)))
# cent_2050_Median <- as.data.frame(st_centroid(SDMshapes_2050_Median))
# cent_2050_Median<-cbind(cent_2050_Median,st_coordinates(st_centroid(SDMshapes_2050_Median)))
# cent_2090_Median <- as.data.frame(st_centroid(SDMshapes_2090_Median))
# cent_2090_Median<-cbind(cent_2090_Median,st_coordinates(st_centroid(SDMshapes_2090_Median)))
# 
# # cent_2050_Q10 <- st_centroid(mypolies_2050_Q10)
# # cent_2090_Q10 <- st_centroid(mypolies_2090_Q10)
# # cent_2050_Q90 <- st_centroid(mypolies_2050_Q90)
# # # cent_2090_Q90 <- st_centroid(mypolies_2090_Q90)
# 
# for (sp in spp){
# 
#   n<-which(spp==sp)
#   nc<-which(gsub(" ","_",cent_Current$SPECIES)==sp)
#   n2050<-which(gsub(" ","_",cent_2050_Median$SPECIES)==sp)
#   n2090<-which(gsub(" ","_",cent_2090_Median$SPECIES)==sp)
#   
#   DatSumMod$lat_Current[n] <- cent_Current$Y[nc]
#   DatSumMod$lat_2050[n] <- cent_Current$Y[n2050]
#   DatSumMod$lat_2090[n] <- cent_Current$Y[n2090]
#   
#   DatSumMod$long_Current[n] <- cent_Current$X[nc]
#   DatSumMod$long_2050[n] <- cent_Current$X[n2050]
#   DatSumMod$long_2090[n] <- cent_Current$X[n2090]
#   
# }
# 
# #set up variables:
# R<- 6371e3
# latpi1 = DatSumMod$lat_Current * pi/180
# latpi2 = DatSumMod$lat_2050 * pi/180
# latpi3 = DatSumMod$lat_2090 * pi/180
# longpi1 <- long_Current * pi/180
# longpi2 <- long_2050 * pi/180
# longpi3 <- long_2090 * pi/180
# 
# Dlatpi2050 = (lat2-lat1) * pi/180;
# Dlatpi2090 = (lat3-lat1) * pi/180;
# Dlongpi2050 = (long2-long1) * pi/180
# Dlongpi2090 = (long3-long1) * pi/180
# 
# #calculate distances:
# a2050 <- sin(Dlatpi2050/2) * sin(Dlatpi2050/2) +  cos(latpi1) * cos(latpi2) *  sin(Dlongpi2050/2) * sin(Dlongpi2050/2)
# a2090 <- sin(Dlatpi2090/2) * sin(Dlatpi2090/2) +  cos(latpi1) * cos(latpi3) *  sin(Dlongpi2090/2) * sin(Dlongpi2090/2)
# c2050 <- 2 * atan2(sqrt(a2050), sqrt(1-a2050))
# c2090 <- 2 * atan2(sqrt(a2090), sqrt(1-a2090))
# DatSumMod$DistM_2050 <- R * c2050 #distance in meters between current and 2050 median distribution
# DatSumMod$DistM_2090 <- R * c2090 #distance in meters between current and 2090 median distribution
# 
# #calculate directions:
# DirRad2050 = atan2(sin(Dlongpi2050) * cos(latpi2) , cos(latpi1)*sin(latpi2) - sin(latpi1)*cos(latpi2)*cos(Dlongpi2050))
# DirRad2090 = atan2(sin(Dlongpi2090) * cos(latpi3) , cos(latpi1)*sin(latpi3) - sin(latpi1)*cos(latpi3)*cos(Dlongpi2090))
# DatSumMod$DirDeg2050 = DirRad2050*180/pi #direction in degrees
# DatSumMod$DirDeg2090 = DirRad2090*180/pi #direction in degrees
# # + = east, - = west
# # 0 = north, 90 = east, -90 = west, 180/-180 =south
# 

# 
# #distance on sphere between point 1 & 2
# 
# R<- 6371e3
# latpi1 = lat1 * pi/180
# latpi2 = lat2 * pi/180
# longpi1 <- long1 * pi/180
# longpi2 <- long2 * pi/180
# 
# Dlatpi = (lat2-lat1) * pi/180;
# Dlongpi = (long2-long1) * pi/180
# 
# a <- sin(Dlatpi/2) * sin(Dlatpi/2) +  cos(latpi1) * cos(latpi2) *  sin(Dlongpi/2) * sin(Dlongpi/2)
# c <- 2 * atan2(sqrt(a), sqrt(1-a))
# DistM <- R * c #distance in meters
# 
# #direction on sphere between point 1 & 2
# 
# DirRad = atan2(sin(Dlongpi) * cos(latpi2) , cos(latpi1)*sin(latpi2) - sin(latpi1)*cos(latpi2)*cos(Dlongpi))
# DirDeg = DirRad*180/pi #direction in degrees
# # + = east, - = west
# # 0 = north, 90 = east, -90 = west, 180/-180 =south

